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SUMMARY 


The  project  is  aimed  at  seismic  surveillance  as  part  of  on-going  efforts  for  improving 
nuclear  test  ban  verification  capabilities.  The  problem  is  complex  in  the  sense  that 
underground  explosions  are  most  efficiently  monitored  by  seismic  means,  but  that  the 
distinction  between  signals  emitted  by  natural  earthquakes  and  explosions  remains  unclear. 
In  other  words,  seismic  wave  propagation  in  a  heterogeneous  Earth  may  easily  mask 
specific  source  signatures. 

In  section  2  we  present  new  techniques  for  fast  and  robust  event  location.  In  case  of 
available  arrival  times  from  a  network  of  seismograph  network  (aperture  up  to  10-15°)  a 
variant  of  Geiger’s  method  for  fitting  arrival  times  to  travel  time  tables,  gave  average 
location  errors  less  than  1°  in  the  teleseismic  distance  range.  In  another  study,  the  location 
schemes  are  tied  to  the  slowness  vector  as  easily  derived  from  arrays  and  3-component 
stations,  namely  azimuth  error  minimization  and  slowness  vector  summation  on  a  sphere 
for  N  arbitrary  positioned  stations.  These  new  techniques  have  already  proved  useful  for 
fast  event  location,  and  thus  help  to  sort  out  the  numerous  reportings  from  the  many 
seismograph  stations  distributed  globally. 

In  section  3  we  deal  with  mapping  of  crustal  structures  in  various  parts  of  Fennoscandia 
where  the  3  seismic  arrays  NORESS,  ARCESS  and  FINESA  are  sited.  The  starting  point 
is  the  NORESS  site  area  using  relative  Pn-arrival  times  in  combination  with  a  modified 
ACH  inversion  scheme.  Due  to  the  small  aperture  of  NORESS  the  mapping  depth  range 
is  also  of  the  same  order,  namely  3  km.  The  major  result  is  that  the  area  is  relatively 
homogeneous  (rhyolite/granite)  except  for  an  intensive  gabbro  body  in  the  NE  quadrant 
which  appears  to  be  roughly  mushroom  shaped  thinning  towards  the  array  center.  In  the 
next  study,  we  have  used  all  available  seismic  reflection  and  refraction  profiling  results  for 
S.  Scandia  so  as  to  produce  a  crustal  thickness  map  for  this  area.  The  crustal  thinning 
associated  with  the  Oslo  Rift  is  pronounced,  of  the  order  of  3-5  km,  and  thus  may  cause 
some  of  the  regional  wavefield  complexities  as  observed  at  the  nearby  NORESS  array. 

Observed  arrival  times  of  local  seismic  crustal  phases  like  Pg,  Sg  (Lg)  and  sub-crustal 
chases  Pn  and  Sn  from  the  Fennoscandian  seismograph  network  are  easily  extracted  from 
computer  files  at  the  national  seismograph  center  in  Helsinki.  Using  a  conjugate  gradient 
scheme  we  have  made  a  tomographic  estimate  of  sub-Moho  lateral  velocity  variation  along 
with  event  hypocenter  parameters  presuming  crustal  thicknesses  to  be  known.  The  most 
prominent  P  and  S  velocity  anomalies  are  found  in  the  northwest  (in-land  Lofoten)  and 
southwest  (M0re)  of  Norway,  and  for  shear  waves  also  in  central  Finland.  Structural 
information  of  the  above  kind  is  essential  for  wavefield  modelling  of  the  mentioned  array 
recordings.  For  example,  some  of  the  observed  anomalies  in  array  azimuths  and  apparent 
velocities  might  be  explained  by  structural  information  as  opposed  to  purely  empirical 
corrections.  Anyway,  the  results  reported  here  arc  to  our  knowledge  the  most  detailed 
published  so  far  for  this  region. 

In  section  4  NORESS  P-wave  coda  of  75  s  duration  from  8  teleseismic  events  of  widely 
different  azimuths  have  been  examined  using  both  array  and  3C  analysis  techniques. 
Although  source-end  scattering  could  not  be  separated  from  the  source  pulse  per  se,  we 
favor  the  hypothesis  of  long  source  duration  supported  by  observed  slow  beam  amplitude 
decay  rates.  The  majority  of  receiver-end  scattering  contributions  appears  to  be  P-to-Rg 
conversions  in  both  forward  and  backward  modes  from  two  nearby  areas  with  pronounced 
topographic  reliefs,  namely  Bronkeberget  (Dist~10  km,  Azi~80°)  and  Skreikampen  (Dist~30 
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km,  Azi~225°).  The  scattering  is  multiple  in  the  sense  that  both  the  primary  and  secondary 
phases  from  the  source  region  contribute  to  the  Rg-scattering.  P-to-S  scattering  constitute 
a  significant  part  of  the  receiver-end  scattering  in  the  later  part  of  the  coda  but  is  more 
diffuse  than  Rg.  P-to-P  scattering  is  weak  and  mostly  confined  to  the  immediate  vicinity 
of  the  array.  These  results  are  interesting  on  two  accounts,  firstly  we  have  pin-pointed  two 
specific  scattering  sources  contributing  significantly  to  the  P-wave  coda  at  the  NORESS 
array  and  secondly  by  invoking  the  principle  of  reciprocity  we  may  expect  similar 
scattering  mechanisms  to  be  responsible  for  scattering  near  the  seismic  source. 


2 


Geophys.  J.  hit.  (1990)  100,  No.  3 


Teleseismic  epicentre  locations  from  arrival  times  at  regional 
networks 

B.  O.  Ruud 

Geology  Department,  University  of  Oslo,  PO  Box  1047  BUndern,  N-0316  Oslo  3.  Norway 


SUMMARY 

Arrival  time  reportings  from  local  and  regional  networks  can  be  used  for 
preliminary  epicentre  location  of  teleseismic  events  via  slowness  and  azimuth 
estimation.  Robust  techniques  for  doing  so  are  presented,  and  their  relative 
performance  tested  on  ISC  listed  P  arrival  times  of  Fennoscandian  stations  (aperture 
ca.  10°).  For  the  10  events  analysed,  the  most  efficient  technique,  a  variant  of 
Geiger’s  method  for  fitting  arrival  times  to  traveltime  tables,  gave  average  location 
errors  of  less  than  1°  for  distances  from  15°  to  85°. 

Key  words:  arrival  times,  epicentres,  teleseismics. 


% 


» 


1  INTRODUCTION 

The  increasing  daily  flow  of  arrival  times  readings  into 
seismological  centres  complicates  the  process  of  pregrouping 
into  /‘vpni  families  necessary  for  nr'*cise  earthq"’Ve 
location.  This  work  is  much  facilitated  when  preliminary 
epi  entre  locations  are  available  from  national  and  local 
cen  -s.  While  epicentre  locations  are  often  reported  for 
local  events,  distant  phases  are  usually  given  with  no 
indication  of  source  location.  The  phase  association  problem 
for  globally  distributed  stations  can  thus  be  very  difficult, 
but  for  a  regional  subset  of  stations  the  problem  is  easier  to 
handle  because  arrival  times  for  a  given  event  will  be  within 
a  shorter  time  interval.  This  requires  the  solution  of  a 
well-known  problem  in  network  and  array  operation,  i.e. 
that  of  locating  earthquakes  based  on  azimuth  and  slowness 
estimates  (see  also  Cassidy  el  al.  1989).  In  this  papjer,  we 
investigate  methods  capable  of  giving  robust  preliminary 
epicentre  locations  from  arrival  times  of  stations  distributed 
over  an  area  of  aperture  of  about  10°  or  less. 

2  METHODS  FOR  NETWORK  SLOWNESS 
VECTOR  ESTIMATION 

The  problem  of  locating  earthquakes  far  outside  the  area  of 
available  stations  is  generally  very  unstable  using- only  P 
phases  as  demonstrated  by  Buland  (1976)  and  Thurber 
(1985).  This  is  due  to  a  lack  of  data  that  makes  the  normal 
equations  almost  singular.  However,  by  assuming  a  fixed 
depth  (or  by  determining  source  depth  by  depth  phases)  the 
problem  is  managable  in  the  range  15°-85°  because 
distances  in  this  range  have  a  one-to-one  correspondance 
with  slowness  (see  Fig.  la).  For  other  distances,  the  first 


arrivals  will  be  P„,  P^-^ft,  or  PKP,  all  with  almost  constant 
slowness  which  ensures  poor  distance  resolution. 
Differential  traveltimes  like  PP-P  or  S-P  can  give  precise 
distance  e.stimates  (Ruud  et  al.  1988),  but  such  methods  are 
not  discussed  here. 

The  idea  of  using  azimuth  and  slowness  estimates  from 
arrival  times  at  a  local  seismic  network  is  not  new,  see  for 
example  Kelly  (1964)  and  Otsuka  (1966).  The  simplest 
approach  is  to  assume  a  plane  wavefront  and  a  local 
flat-earth  approximation  and  solve  for  the  least-squares 
estimate  of  the  horizontal  slowness  vector. 

2.1  Method  1:  plane  wavefront 

The  equation  for  a  plane  wavefront  as  a  function  of  latitude 
(p  and  longitude  A  can  be  written 

t(<P,  A)  =  t{(pQ,  Ay)  +  (A  —  A,|)  cos  tpvU^  -(-  {ip  —  (p^f)uf,|  (1) 

where  and  the  east  and  north  components  of  the 

horizontal  slowness  vector  and  (<p„,  Aq)  is  a  reference  pwint 
in  the  array.  Since  the  equation  is  linear  in  the  unknowns 
i(,<Po>  An),  Ue  3nd  Uf^,  standard  linear  inverse  methods  can 
be  applied.  Azimuth 

y„  =  arctan(uE/MN) 

and  slowness 


relative  to  the  reference  point  are  then  found. 

This  approach  is  satisfactory  for  arrays/networks  with 
aperture  up  to  about  100  km,  depending  of  the  accuracy  of 
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1.  (a)  f'-wave  Iraveltir.ie  (solid  curve)  and  slowness  (dashed  curve)  for  surface  focus.  Data  are  from  Herrin’s  (l‘)6)i)  tables,  (b) 
Expected  location  error  due  to  an  uncertainty  in  the  horizontal  slowness  vector  of  0.1  s  degree"'  for  P  phase.  The  lower  curve  is  for  azimuth 
and  the  middle  curve  for  slowness  error.  Hie  upper  curve  is  the  total  location  error. 


the  arrival  time  measurements.  For  larger  arrays,  slowness 
variations  over  the  area  of  the  network,  wavefront  curvature 
due  to  finite  epicentral  distance,  and  distortion  from  the 
fiat-earth  transformation  require  other  schemes. 


2,2  .Method  2:  curved  wavefront  based  on  Taylor 
expansion 

N'arious  methods  have  been  proposed  to  remedy  the 
dcficienccs  of  the  plane  wavefront  technique.  Most  of  these 
are  based  on  repicsentations  of  the  wavefront  as  second  or 
higher  order  polynomial  surfaces  in  Cartesian  or  spherical 
coordinates.  Also  non-integer  powers  of  latitude-ioncitude 
coordinates  have  been  tried  (Husebye  1969). 

Although  these  methods  are  capable  of  representing  many 
kinds  of  curved  surfaces,  the  estimated  surface  will  not 
ncressarily  resemble  the  wavefront  expected  from  a  point 
source  in  a  spherical  symmetric  earth  model.  In  order  to 
obtain  reliiable  estimates,  the  number  of  free  parameters 
should  also  be  kept  as  low  as  possible. 

The  least  ambigious  assumption  to  make  about  the  arrival 
time  IS  that  it  IS  a  function  of  epicentral  distance  only  (for 
fixed  depth).  To  account  for  smooth  variations  in  the 
slowness  over  the  network,  the  wavefront  could  be 
approximated  by  a  second-order  Taylor  expansion  formula 
in  the  arrival  time  /(A): 


f(4)  =  f(A,)  +  ^|  (4 


>  1 


(4 -A.) 


(2) 


In  order  to  calculate  the  distance  for  each  station  a  'trial' 
epicentre,  defined  by  its  distance  Ai  and  azimuth  y„  relative 
to  the  reference  point,  is  used.  Distance  4  is  a  non-linear 
fuivion  of  the  epicentre  coordinates  (Ai,  Tu)  so  an  iterative 
tec:  nique  is  necessary  to  solve  for  a  least-squares  estimate 
of  the  unknowns.  The  epicentre  coordinates  are  updated  for 
each  iteration  and  provided  that  the  algorithm  converges  to 
ti.c  correct  solution,  tiic  method  is  exact  in  azimuth  while  a 
S'cond-order  approximation  is  made  for  distance.  In 
addition  to  the  epicentre  coordinates  there  are  three 


unknown  parameters;  /(A,,),  5r/34|^  ,  and 
Some  of  these  parameters  are  unnecessary  under  certain 
assumptions.  For  instance,  fixed  values  might  be  used  for 
(almost  constant  for  many  phases  and  distance 
ranges)  or  for  Ai  'f  distance  is  large  compared  to  aperture. 

Note  that,  so  far.  no  assumptions  have  been  made  about 
the  typte  of  phase  used  (except  that  it  must  be  the  same  for 
all  stations)  or  the  velocities  in  the  earth  model.  The  next 
step  is  to  use  the  estimated  slowness  and  convert  it  into 
distance  using  traveltime  tables  or  earth  velocity  models. 
Thus,  proper  phase  identification  must  be  made,  but  in 
principle  the  slowness  and  azimuth  estimates  discussed 
above  are  independent  of  phase  type.  Relatively  large 
discrepancies  are  sometimes  found  between  the  parameter 
A»,  depending  on  wavefront  curvature,  and  the  distance 
corresponding  to  the  estimated  slowness.  This  is  not 
unexpected  when  distance  is  large  compared  to  aperture  or 
if  only  few  readings  with  large  errors  are  available.  Since 
distance  found  from  slowness  conversion  is  far  more 
reliable,  4o  is  here  considered  as  a  parameter  necessary  to 
represent  a  realistic  wavefront  and  not  as  the  final  distance 
cstim’.te. 


2.3  Method  3:  fitting  arrival  times  directly  to  traveltime 
tables  or  earth  models 

In  many  cases,  the  first  P  phase  is  easily  identified  as  such 
and  the  possible  advantage  of  a  phase-indepiendent  slowness 
estimate  is  not  relevant.  Also,  the  known  coupling  between 
4,  3f/54  and  ci^tldd}  due  to  the  information  accumulated 
in  traveltime  tables  (or  earth  velocity  models),  should  be 
utilized  to  reduce  the  number  of  unknowns.  This  could  be 
done  by  adding  certain  constraints  to  the  model  and  solving 
for  just  one  of  these  three  distance-depiendent  parameters. 
If  we  choose  to  use  r(4„).  A,,  and  y„  as  unknowns,  we  arc  in 
fact  back  to  the  traditional  earthquake  location  scheme 
originally  due  to  Geiger  (1910).  The  equation  for  the 
wavefront  can  then  be  written  as 

t(4)  =  /o+r(4)  (3) 


# 
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wliere  r,,  is  or  omc  and  7(4)  denotes  the  travciiime 
derived  fror  •  les  or  an  earth  model.  As  with  method  2, 
distance  jepends  on  the  epicentre  coordinates  given  by 
distance  4,,  and  azimuth  y„  relative  to  the  reference  point. 
Tnis  representation  is  preferred  because  it  seems  to  improve 
convergence  compared  to  sphencal  coordinates  when 
sohini’,  this  non  linear  inverse  problem  with  a  Gauss- 
Newton  Iteration  method  (see  also  Thurber  1985). 
Appendix  A  shows  how  the  partial  derivatives  of  f  arc 
calcula'ed  with  this  epicentre  representation. 


2.4  Error  analysis 

Formal  error  analysis  can  be  carried  out  as  usual  for 
least-squares  estimation.  What  is  given  here  is  only  a 
simplified  analysis,  with  emphasis  on  how  errors  propagate 
from  observations  via  azimuth  and  slowness  into  location 
errors  Discussion  of  array  location  capabilities  can  also  be 
found  in  Gjoystdal,  ifusebye  &  Rieber-Mohn  (1973). 

1  he  accuracy  of  the  estimated  slowness  vector  depends  on 
the  array  geometry  and  the  accuracy  of  the  arrival  times.  As 
a  rough  approximation,  the  slowness  uncertainty  &u  is  equal 
to  the  time  uncertainty  divided  by  the  aperture,  (i.e.  the 
diameter  of  the  least  circle  surrounding  all  stations) 
provided  that  the  network  is  not  too  elongated.  For 
instance,  a  network  of  10°  aperture  with  relative  arrival  time 
errors  of  1  s,  would  give  slowness  errors  of  the  order  of  0.1  s 
decree"'.  Precision  may  improve  with  many  stations,  but 
aperture  is  still  the  main  parameter  controlling  slowness 
accuracy 

.Azimuth  errors  tire  slowness  dependent  and  is  greatest  for 
near  vertical  incidence  (low  slowness)  The  azimuth 
uncertainty  by  is  given  by 


(in  radians).  Thus,  for  a  certain  phase,  the  expected  azimuth 
error  is  a  function  of  distance  through  the  slowness.  The 
location  error  due  to  an  azimuth  error  is  scaled  by  a 
distance  dependent  term; 

=  sin  4  6)'. 

When  slowness  is  converted  into  distance,  the  distance 
uncertainty  is  controlled  by  the  phase-dependent 

parameter  d'UdA^,  i.e.  the  slope  of  the  slowness  curve  in 
Fig  1(a).  The  relation  here  is 

he  .=  — ; - v  . 

|<9-//a41 

The  total  epicentre  location  error  will  then  be 
''’<’101  *  \43ej  -t-  bCy . 

In  Fig.  Kb),  de^,  and  Ae,,,,  are  shown  for  the  P  phase. 
Computations  were  made  using  Herrin  tables  and  =  0. 1  s 

degree  '. 

In  addition,  there  are  bias  errors  due  to  velocity 
anomalies.  This  is  a  whole  subject  on  its  own  and  may  often 
be  the  final  factor  controlling  location  precision.  However, 
the  effect  is  likely  to  be  most  serious  for  small  arrays,  and 
can  to  a  large  degree  be  corrected  for  by  introducing  station 
delays  and/or  regional  location  callibrations.  For  instance. 


the  average  slowness  vector  deviation  at  the  NORSAR 
array  was  reduced  from  (I.?  to  (Ms  degree"'  using  an 
interpolation  scheme  based  on  93  'master  events'  (Berteus- 
sen  1976). 

3  RE:SULTS  from  real  data  TEST.S 

The  three  methods  have  been  tested  on  arrival  times  listed 
in  ISC  bulletins.  A  total  of  10  events  within  appropriate 
distance  ranges  and  with  a  reasonable  number  of  /’-phase 
readings  for  stations  in  the  Fennoscandia  were  selected 
randomly.  The  epicentre  solutions  as  reported  by  ISC  are 
listed  in  Table  1,  and  the  stations  used  for  each  event  are 
given  in  Table  2.  Testing  was  first  attempted  using  all 
reporting  stations  in  Norway,  Sweden,  Finland  and 
Denmark,  but  after  some  experiments  it  was  (ound  that  the 
results  actually  improved  by  restricting  slaiioris  only  to 
Sweden  and  Finland,  constituting  the  core  of  the  Baltic 
shield,  which  is  likely  to  be  a  more  laterally  homogeneous 
region  (see  also  Noponen  1974)  The  stations  are  distnbuted 
in  the  ranges  60°-7U‘’N  and  12°-2S°Li.. 

The  plane  wave  solution  was  used  as  a  starting  point  tor 
the  Iterative  non  linear  least-squares  estimation  of  method  2 
and  3.  Weighting  and  damping  according  to  reasonable  a 
priori  information  as  described  by  Jordan  &  Sverdrup  (19S1) 
and  Jackson  &  Matsu'ura  (19S.S)  were  used  to  stabilize  the 
inversion.  Convergence  was  assured  through  the  use  of 
so-called  spacer  gradient  steps  (see  Luenberger  1973). 

For  each  event,  the  reference  point  was  taken  as  the 
arithmetic  mean  of  the  station  coordinates,  and  the 


Table  1.  Events  used  m  lesls.  J  tie  regions  of  the  events  are;  ! 
Greece,  2.  Carlslierc  Ridge,  Kodiak  Islands,  4.  Pakistan,  5 
Kushu.  Japan.  6.  Mananas  Islands.  7;  'radzhikisian,  USSR.  8 
Greece.  9;  Afghanistan,  10.  Equador. 
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Table  2.  Stations  used  in  analysis.  These 
arc  all  stations  in  Sweden  and  Finland 
listed  in  the  bulletins  for  the  selected 
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Table  3.  Estimated  azimuths,  slownesses  and  distances  for  the 
Table  1  events.  The  (•  line  for  each  event  gives  the  azimuth 
and  distance  of  the  epicentre  relative  to  the  reference  point.  The 
slowness  is  found  from  Herrin’s  (1968)  tables  for  given  distance 
and  source  depth.  For  method  1  and  2,  the  slowness  is  estimated 
directly  and  the  corresponding  distances  are  found  from  tables 
assuming  surface  focus,  while  for  method  3  it  is  (he  other  way 
round  err  ’  is  the  total  epicentre  location  error  using  the 

estimated  azimuths  and  distances,  and  ‘RMS'  is  (he  root  mean 
square  of  the  time  residuals  after  estimation.  For  the  agencies 
rms  residuals  have  been  calculated  from  the  bulletin  residuals 
after  subtracting  an  average  for  each  event. 
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estimated  azimuths,  slownesses  and  distances  relative  to  this 
point  are  given  in  Table  3. 

As  expected,  it  is  clear  from  the  large  location  errors  and 
rms  residuals  of  method  1  that  the  plane  wave  and  flat-earth 
approximation  is  not  suitable  for  a  network  of  this  size.  Still, 
azimuth  estimates  are  quite  good  (errors  less  than  3°),  while 
slowness  errors  seems  to  be  azimuth  dependent  (small  to  the 
north  and  south  and  larger  to  the  east  and  west). 

For  most  events,  method  2  gives  very  good  estimates  with 
small  rms  residuals.  The  largest  location  error  is  for  event  4 
which  also  has  the  largest  rms.  In  this  case,  the  parameters 


d^i/dA~\^^  and  attain  unlikely  values  when  adjusting  to 
the  apparently  poor  observational  data.  Tlius,  method  2 
could  be  unstable  for  data  with  large  errors. 

As  expected,  method  3  handles  the  'erroneous'  data  of 
event  4  better  than  method  2  due  to  fewer  free  parameters. 
Also  this  method  accounts  for  variations  in  d~t/dA~  over  the 
network,  most  important  for  distances  less  than  30°.  The 
largest  location  error  here  is  for  event  10  and  is  due  to 
distance  error.  Fig.  1(b)  indicates  that  this  is  to  be  expected 
for  distances  beyond  85°.  The  location  error  for  event  9  is 
mainly  due  to  the  fact  that  this  is  an  intermediate  depth 
earthquake.  As  seen  from  the  table,  the  slowness  estimate 
here  is  actually  better  than  for  method  2,  although  the  latter 
coincidentally  gives  a  better  epicentre  location  when  zero 
depth  is  assumed.  The  largest  errors  in  slowness  are  found 
for  distances  less  than  30°  and  could  be  due  to  regional 
variations  in  the  upper  mantle  P  velocities.  The  Herrin 
tables  u.sed  are  adjusted  to  fit  traveltimes  of  the  central  US 
for  regional  distances  (Herrin  1968). 

4  CONCLUSIONS 

Using  F  arrival  times  reported  from  a  network  w'ith  aperture 
less  than  10°  we  are  able  to  locate  epicentres  to  within  1.5° 
in  the  distance  range  15°-85°.  Possible  exceptions  are  deep 
earthquakes  or  distances  close  to  30°  (sec  Fig.  lb). 
Estimates  can  probably  be  improved  by  using  empirical  time 
corrections  or  location  callibrations  not  attempted  here. 

The  largest  contribution  to  the  location  error  comes  from 
the  distance  estimate.  As  azimuth  is  expected  to  be  very 
accurate  (within  1°)  for  all  distances  and  phases,  except  core 
phases,  other  information  could  be  utilized  to  resolve 
distance  ambiguity.  This  includes  seismicity  information  and 
difl'erential  traveltimes. 

Of  the  three  methods  tested,  a  variant  of  Geiger’s  method 
for  fitting  arrival  times  directly  to  traveltime  tables  gave  the 
most  stable  results.  This  method  relies  on  proper  phase 
identification,  so  the  phase-independent  method  2  may  be 
valuable  for  initial  investigations.  The  plane  wave  method  is 
obviously  only  suited  for  very  small  networks/arrays,  but 
may  in  this  case  be  preferable  due  to  its  simplicity. 
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APPENDIX 

In  order  to  apply  the  Gauss- Newton  iteration  technique  for 
the  least-squares  estimation  in  method  2  and  3,  it  is 
necessary  to  calculate  the  partial  derivatives  of  i{A)  with 
respect  to  the  parameters  to  be  estimated.  As  distance  A  is  a 
function  of  A„  and  y„,  we  have  (for  method  3) 

di  dT  <9 A 
dA(y  dA  5Ao 
and 

Bt  dT  BA 

Denoting  the  distance  and  azimuth  from  'he  reference  point 
to  the  /'th  station  as  d,  and  z,  resfjectively,  we  find  by 
considering  the  triangle  formed  by  the  epicentre,  the 
reference  point  and  this  station 

cos  A  =  cos  Ao  cos  d,  -I-  sin  A^  sin  d,-  cos  (z,  -  Vq)-  (A2) 

Deriving  this  with  respect  to  A<,  and  y„  we  find 
3  A 

- —  =  (sin  A, I  cos  d,  —  cos  Ay  sin  d,  cos  (z,-  -  y(|)]/sin  A  (A3) 
dA„ 

and 

aA 

—  =  -sin  A„  sin  d,  sin  (z,  -  yo)/sin  A.  (A4) 

SYo 

These  relations  are  also  used  when  calculating  the 
derivatives  for  method  2. 
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ROBUST  AND  RELIABLE  TECHNIQUES  FOR  EPICENTER  LOCATION 
USING  TIME  AND  SLOWNESS  OBSERVATIONS 

By  Frances  Cassidy,  Anders  Christoffersson,  Fystein  S.  Husebye, 

AND  Bent  O.  Ruud 

Abstract 

The  ever-increasing  flow  of  parameterized  and  waveform  data  into  various 
kinds  of  seismological  centers  cannot  be  managed  property  unless  preliminary 
epicenter  locations  are  available.  Here  we  demonstrate  robust  and  flexible 
techniques  for  fast  and  reliable  event  location  using  the  P  slowness  vector,  which 
is  easily  derived  from  arrival  times  and/or  waveform  data  from  arrays,  networks, 
and/or  single-site  three-compo  ient  stations.  Location  schemes  are  tied  to  1) 
azimuth  minimization  and  2)  slowness  vector  summation  on  a  sphere  using 
slowness  from  N  arbitrarily  positioned  stations.  The  advantage  of  using  azimuth 
alone  is  that  no  assumption  is  needed  of  phase  type,  distance  range,  or  structure 
(travel-time  tables).  The  viability  of  our  location  techniques  are  demonstrated 
using  a  variety  of  P  recordings  from  networks  and  three-component  stations. 

Introduction 

The  problem  of  precisely  locating  earthquakes  has  intrigued  seismologists  for 
decades  and  still  does.  The  various  approaches  published  are,  with  few  exceptions, 
related  to  Geiger’s  (1910)  iterative,  least-squares  method  on  the  basis  of  P-wave 
travel-time  obser\'ations  (e.g.,  see  Jordan  and  Sverdrup,  1981;  Pavlis,  1986).  Notable 
extensions  of  the  Geiger  scheme  are  those  of  Julian  (1973)  who  also  included 
slowness  estimates  from  array  stations,  and  those  of  Aki  and  Lee  (1976)  and 
Spencer  and  Gubbins  (1980)  who  devised  novel  techniques  for  jointly  estimating 
hypocenter  and  earth  structure  parameters.  Arrays  and  single  three-component 
stations  can  provide  fast  although  less  precise  epicenter/hypocenter  locations  (e.g., 
see  Bratt  and  Bache,  1988;  Ruud  et  ai,  1988). 

The  topic  of  this  paper  is  that  of  novel  schemes  for  automated,  preliminary  event 
locations  on  the  basis  of  P-wave  slowness  vector  estimates  from  local/regional 
networks  of  conventional  stations,  arrays,  and/or  digital  three-component  stations. 
It  reflects  a  rather  pressing  problem  at  many  seismological  centers,  namely,  an 
increasing  daily  flow  of  P  arrival  reportings  from  regionally  and  globally  distributed 
stations  from  which  grouping  into  event  families  is  a  prerequisite  to  undertaking 
formal  event  location  procedures.  This  process  is  by  no  means  trivial  and,  in 
addition,  is  work-intensive  but  is  much  facilitated  for  events  for  which  preliminary 
epicenter  locations  are  available.  An  illustrative  example  here  is  NEIC/USGS  which 
receives  daily  between  1000  to  2000  P  arrival  readings,  out  of  which  only  about  60 
per  cent  are  associated  with  specific  hypocenter  solutions.  The  relatively  low 
association  percentage  partly  reflects  a  formal  constraint  of  four  reportings  for 
defining  an  event  location.  Today,  an  increasing  number  of  ‘smart’,  digital  three- 
component  stations  are  deployed  on  a  global  scale.  By  ‘smart’  we  mean  that  real¬ 
time  processing  is  used  for  the  detection  of  incoming  signals  and  the  novel  three- 
component  processing  schemes  of  Christoffersson  et  al.  (1988)  and  Roberts  et  al. 
(1989)  could  easily  be  incorporated  for  extracting  the  P  slowness  vector  as  well. 
Note  that  the  epicenter  location  techniques  developed  here  are  tested  on  signal 
parameters  derived  off-line,  as  demonstrated  in  Figure  1. 
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Fig.  1.  Extraction  of  Slowness  Vector.  Signal  from  the  magnitude  m  =  6.0  Greek  earthquake  (Event 
3,  Table  2),  recorded  at  the  3-component  station  at  Chaingmai  (18.8N,  99.0E).  The  plotted  contours  are 
probability  levels  as  a  function  of  time  and  azimuth  using  a  0.50  min  time  window  (WLS)  and  a  time 
step  (DS)  of  0.25  min  The  data  are  velocity  filtered  (CUTL.  CUTU)  from  5-55  km/sec,  and  the 
probability  threshold  in  filtering  the  original  traces  (SYNT)  is  0.5.  TIME,  AZIM,  PROB,  and  VELO 
values  read  at  cross  in  figure  are  given  at  the  right.  Also  shown  are  estimated  latitude  and  longitude 
(LAT,  LONG)  and  travel  time  (TRAV).  Other  notations;  AV  =  smoothing  parameter;  AL  =  azimuth 
interval  in  degrees;  SRATE  =  sampling  rate  (min). 


Data  Availability 

Before  presenting  and  discussing  techniques  for  fast  epicenter  locations,  a  few 
words  on  the  data  at  hand  for  this  task  are  necessary.  We  differentiate  here  between 
observations  provided  by  local/regional  networks,  arrays,  and  three-component 
stations,  although  the  common  denominator  here  is  estimates  of  the  P-wave 
slowness  vector,  S.  Note  that  in  practice  S  is  replaced  by  its  corresponding  epicentral 
distance  via  standard  travel-time  tables.  Furthermore  the  S  vector  often  does  not 
have  a  stringent  definition  and  mav  be  interpreted  as  apparent  velocity,  dT/dl  or 
distance. 

1.  Local/regional  networks  (aperture  ]  to  15  degrees):  Reported  P-wave  param¬ 
eters  are  arrival  time  T,,  amplitude  and  period  for  the  ith  station.  T,  observa¬ 
tions  are  transformed  to  an  S  estimate  using  a  modified  Husebye  (1969) 
scheme  as  detailed  by  Ruud  (1990). 
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2.  Arrays:  Same  as  above  but  the  S  estimate  is  either  tied  to  the  detecting  beam 
location  or  obtained  directly  from  f-k  analysis  or  variants  thereof. 

3.  Three-component  stations:  Same  as  above  but  S  is  obtained  via  an  estimate 
of  the  P-wave  incidence  (polarization)  angle. 

The  precision  in  S-estimat.es  depends  on  network/array  aperture  and,  for  three- 
component  stations,  on  wavelength.  For  small  arrays  apparent  wavelength  is  im¬ 
portant  and  should  preferably  be  less  than  aperture  (Harris,  1982).  For  local, 
regional,  and  core  phases  the  epicentral  distance  estimates  corresponding  to  S 
conversions  via  travel-time  tables  are  often  poor  due  to  small  time  gradients  in 
these  ranges.  If  secondary  phases  can  confidently  be  identified,  differential  travel 
times  like  dT(Pn  -  Pg),  dTiPKP  —  PP),  etc.,  would  provide  accurate  distance 
estimates  at  the  above  troublesome  ranges.  For  a  three-component  station,  a  1° 
error  in  vertical  incidence  angle  is  roughly  equivalent  to  a  4°  error  in  distance  at 
teleseismic  ranges.  Azimuth  errors  scale  differently  and  peak  at  90°  distance. 

In  the  above  list  of  observational  data  we  have  included  S  from  digital  three- 
component  stations,  as  this  parameter  can  easily  and  reliably  be  extracted  using  an 
appropriate  analyzing  scheme  (see  Fig.  1).  However,  Harris  (1982)  and  Bratt  and 
Bache  (1988)  consider  such  observations  of  marginal  interest  due  to  a  claimed  poor 
precision  in  S  estimates,  even  in  comparison  to  small  aperture  (3  km)  arrays  like 
NORESS.  It  suffices  here  to  refer  to  a  recent  study  by  Ruud  et  al.  (1988)  who  found 
that  slowness  (azimuth)  estimates  derived  from  NORESS  and  single  three- 
component  station  recctrds  were  of  similar  precision  for  even  small  events 
(SNR  >  2)  at  local/regional  distances.  At  teleseismic  ranges  small  arrays  suffer 
from  increases  in  apparent  wavelength,  while  for  three-component  stations  the 
horizontal  energy  decreases  but  for  adequate  SNRs  do  not  impair  the  S-estimate 
per  se. 


Fast  Epicenter  Location  Techniques 

In  this  section  we  present  event  location  schemes  on  the  basis  of  extracted 
P-wave  signal  parameters  as  outlined  above.  Epicenter  locations  for  azimuth  obser¬ 
vations  alone  are  described  first,  and  then  slowness  vector  location.  However,  before 
doing  so,  we  comment  on  direct  slowness  estimation  for  P  arrival  times  as  reported 
by  local/regional  networks. 

Network  Arrival  Times — Wavefront  Fitting 

Using  spherical  coordinates,  polynomial  wavefront  fitting  is  feasible  for  P  arrival 
times  across  a  network  of  stations  of  aperture  say  in  the  range  of  1  to  15°  and  hence 
by  derivations  can  be  estimated.  The  advantage  of  using  networks  for  S  estimation 
is  that  their  large  apertures  give  precise  results  with  both  azimuth  and  distance 
errors  within  a  few  degrees.  A  potential  problem  is  that  individual  stations  occa¬ 
sionally  report  different  phases,  say  Pg  and  Pn  at  local  distances  and  PKP(DF)  and 
PKP(BC)  in  the  core  window. 

Azimuth  Location — Angle  Minimization 

The  idea  here  is  to  minimize  the  square  differences  between  observed  azimuth 
(a.)  at  N  stations  and  those  tied  to  the  unknown  epicenter  location  {$,  X),  namely; 

N 

F  =  [{a,  -  GAi  ,  X)V.  (1) 

1^1 
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G,  is  defined  as; 


G,  =  sin 


,  sin(^)sin(X,) 
sin(A,) 


(2) 


where  A,  is  the  epicentral  distance  for  the  ith  station. 

Estimates  of  the  standard  error  of  the  epicenter  estimate  can  be  derived  from  its 
covariance  matrix  (a  j  reflects  the  variance  of  the  observed  azimuths); 


cov((^,  X)  ~  “  2) 


(3) 


where  asterisk  denotes  transpose; 


A' 


I 


A. 

A, 

A 

B, 

and 


(4) 


>1. 


di 


,  dG, 
and  B,  =  — ^  . 

d\ 


(5) 


Details  on  dGJdi  and  dGJdX  are  given  in  the  Appendix. 

The  advantage  of  using  the  azimuth  part  of  S  is  simply  that  this  parameter  is 
independent  of  phase  type  and  thus  distance  range  and  travel-time  tables.  In  other 
applications,  azimuth  from  mainly  two  arrays  have  been  used  in  conjunction  with 
arrival-time  observations  in  Geiger-type  inversion  schemes  for  estimating  <t>,  X  (e.g., 
see  Gjdystdal  et  al.  1973;  Bratt  and  Bache,  1988).  Our  preference  for  using  azimuth 
alone  from  several  stations  (three-component  recordings)  is  that  a  precise  prelimi¬ 
nary  location  is  feasible,  which  subsequently  can  be  used  for  phase  identification 
and  association.  The  latter  is  often  a  problem  at  international  seismological  data 
centers. 


Slowness  Vector  Location — Vector  Summing 

In  addition  to  bearing,  S  contains  information  on  epicenter  distance,  that  is,  after 
conversion  to  distance  via  standard  travel-time  tables.  Such  a  conversion  is  in 
practice  only  valid  in  the  teleseismic  window  due  to  small  gradients  in  the  travel- 
time  curve  at  other  distance  ranges  unless  additional  information  is  provided.  Now, 
the  strategy  adapted  for  epicenter  locations  using  slowness  is  by  means  of  vector 
summing  based  on  a  method  developed  by  Fisher  (1953)  for  combining  observations 
of  position  (epicenter)  on  a  sphere  (see  also  Irving,  1964).  The  optimum  event 
location  or  center  of  gravity  for  N  dispersed  epicenter  (tp,  X)  observations  is 
calculated  from  the  vector  sum  of  unit  vectors  of  the  N  observed  directions.  Each 
of  these  directions  is  specified  by  its  3  direction  cosines,  namely; 

North  I  =  cos(0)cos(X); 

South  m  =  cos(<^)sin(X);  (6) 

Vertical  z  =  sin(0). 
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The  best  estimate  of  the  average  direction  is  then  a  vector  sum  of  the  individual 
directions; 


1  N  1  ^  1  ^ 

^  =  P  S  A;  Af  =  p  S  2  =  -  X  2. 

K  ,=  i  n  ,_i  n  ,_i 

where 


2 

'  N 

12 

'  S 

= 

I  4 
.1=1 

+ 

.j“l 

S  2- 

Then,  the  mean  epicenter  location  is  given  by; 

tan(<^)  =  M/L  sin(X)  =  Z. 


(7) 


(8) 


(9) 


Fisher  showed  that  the  probability,  P,  of  a  point  being  located  between  <j)  and 
4>  +  d(i>  is; 


=  o  •  d(t>,  (10) 

2  sin  h{k) 

where  k  determines  the  precision  of  the  points.  The  best  estimate  of  k,  the  so-called 
precision  factor,  k,  is; 


N  -  1 


N- R 


(11) 


Also,  an  accuracy  measure  of  the  mean  direction  can  be  calculated  for  a  particular 
probability  level,  (1  -  P).  The  95  per  cent  confidence  radius  is  given  by; 

140 

and  a  65  per  cent  radius  (corresponding  to  the  standard  error  in  Gaussian  distri¬ 
butions)  is; 


67.5 

These  measures  would  be  useful  in  comparing  locations  based  on  slowness  with 
those  reported  by  NEIC,  besides  serving  as  a  confidence  ellipse  estimate. 

Data  Analysis  and  Results 

In  order  to  demonstrate  the  location  performances  of  the  above  techniques,  S 
vectors  were  extracted  from  three-component  records  at  1  to  90°  distances  using 
the  Christoffersson  et  al.  (1988)  technique.  A  practical  problem  encountered  was 
that  access  to  digital  records  from  well -separated  stations  was  not  easy  and  data 


12 


EPICENTER  LOCATION  WITH  TIME/SLOWNESS  OBSERVATIONS 

quality  was  occasionally  poor.  For  some  stations  the  orientations  of  the  horizontal 
components  appear  to  be  off  by  15  to  20'^.  Nonetheless,  results  to  be  presented 
demonstrate  the  viability  of  our  location  schemes. 

Explosion — Geophone  Record  Analysis 

During  the  Fennolora  seismic  refraction  experiment  (Lund  and  Prodehl,  1985)  a 
number  of  three-component  geophones  were  deployed.  In  addition  to  the  profiling 
shots,  a  nuclear  explosion  at  Semipalatinsk  was  also  recorded  by  chance.  Two  of 
these  events  ha\  e  been  subjected  to  “location”  analysis,  and  details  on  the  outcome 
are  given  in  Table  1. 

Event  1:  Since  epicenter  distances  were  within  a  few  degrees,  only  azimuth 
location  was  feasible.  Despite  an  anomalous  observation  (station  134,  19.4°)  the  net 
mislocation  is  around  6  km,  which  is  considered  very  good  as  the  epicenter  distance 
range  is  200  to  400  km. 

Event  2:  The  aperture  of  the  geophone  line  was  too  small  compared  to  epicentral 
distance  for  basing  the  location  on  azimuth  only.  However,  using  the  slowness 
vectors  the  estimated  epicenter  is  about  1°  off. 

Teleseismic  Events 

We  have  here  analyzed  two  sets  of  observations,  three-component  broadband 
recordings  on  NEIC/ORFEUS  CD  ROMs  and  the  corresponding  P  arrival  times 


TABLE  1 

Fennolora  Profile  Results 


Sution 

Ev<n(  1 

Fennolora  Shm 
(AS47N.  i:„3riKi 

Event  2 

Semipalatinsk  Explosion 
(49-89N.  78,89Ei 

Velocity 

(km/ 

tec) 

Azimuth 

(deni 

dAzi 

(deg) 

Velocity 

(km/ 

see) 

(deg* 

dS 

(degi 

Azimuth 

(deg) 

dAzi 
(deg  < 

Latitude 

I'N’t 

Longitude 

CEi 

14 

7.9 

15.8 

-4.8 

13.5 

42.4 

5.4 

67.3 

-5.4 

494 

88.8 

17 

12.2 

19.8 

-1.9 

13.9 

45.6 

8.7 

79,2 

6.4 

41.1 

84.5 

44 

6.8 

33.6 

2.4 

12.5 

31.0 

-5.6 

75.3 

1.3 

52.4 

70.9 

47 

6.8 

35.6 

2.1 

13.6 

42.6 

6.0 

79.2 

3.8 

43,5 

82.6 

75 

7.1 

83.2 

7.5 

12.9 

35.7 

-0,5 

79.2 

3.8 

48.4 

76.1 

86 

8.4 

110,9 

-0.5 

13.0 

37.6 

1.5 

67.3 

-6.8 

53.0 

85.9 

89 

6.8 

118.9 

0.2 

11.4 

24.6 

-12.6 

67.7 

-9.9 

67.7 

62.1 

95 

7.6 

134.8 

0.7 

13.3 

40.3 

4.2 

65,4 

-11.0 

52.5 

90.9 

110 

7.8 

152.6 

-2.4 

12.8 

34.9 

-1.0 

77.2 

0.2 

50.4 

77.6 

134 

8.7 

148.6 

-19.4 

11.7 

25.6 

-10.0 

71.3 

-6.8 

58.4 

67.9 

146 

8.8 

174.9 

2.6 

signal  clipped 

Estimated  epicenter  Azimuth  (Sec  3.2):  58.43  N  (S.E.  0.04  N);  17.53  E  (S.E.  =  0.12  E);  Slowness 
(Sec  3.3):  50.82  N;  78.85  E  (tv,  =  1.60);  (a,s  =  2.77), 

Three-component  field  geophone  recordings  from  the  Fennolora  deep  seismic  profiling  experiment 
used  in  location  experiments:  Event  1  is  a  Fennolora  shot  (dynamite  charge  of  800  kg),  and  event  2,  an 
underground  nuclear  explosion  at  Semipalatinsk  (USSR)  which  happened  to  be  recorded.  In  event  1. 
only  the  azimuth  part  of  the  S  vector  was  used  for  epicenter  location  as  Pg-  and  Pn -phase  velocities  do 
not  provide  acceptable  distance  estimates.  The  dAzi  column  gives  the  difference  between  the  “true”  and 
observed  azimuths.  For  event  2.  the  S  estimate  for  each  station  provides  an  epicenter  estimate  as  listed, 
but  using  the  vector  summation  scheme  over  ail  reporting  stations,  an  acceptable  epicenter  location 
(error  near  I')  is  obtained.  As  before,  dAzi  and  dA  give  differences  between  “true”  and  obser\’ed  azimuth 
and  distance  estimates,  respectively.  Azimuth  location  is  not  feasible  due  to  the  small  aperture  of  the 
line  array 
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for  Fennoscandian  stations  as  reported  by  ISC.  The  results  in  terms  of  both  azimuth 
and  slowness  locations  are  given  in  Tables  2  and  3  and  for  “networking”  in  Tables 
4  and  5.  The  following  comments  apply:  slowness  observations  from  a  few  three- 
component  stations  suffice  for  making  fast  and  reliable  (within  a  few  degrees) 
epicenter  determinations.  More  refined  estimates  are  feasible  by  introducing  slow¬ 
ness  bias  corrections  and/or  seismicity  information.  Likewise,  networking,  as  dem¬ 
onstrated  here  using  Fennoscandia  P  arrival  times  for  slowness  estimation,  also 
provides  reliable  locations  even  on  a  stand-alone  basis. 


TABLE  2 

Stations  and  Events  used  for  Location  Based  on 
NEIC/ORFEUS  Broadband  Data 


Station 

Event  I 
Ban<U 

Sea 

Event  2 

Event  3 

Event  4 

Ecuaor 

Greece 

Afghaniatan 

ANM 

ANT 

X 

X 

X 

BCA 

X 

CHT 

X 

X 

CTA 

X 

X 

GRF 

X 

X 

GUM 

JAS  ■ 

X 

X 

KON 

X 

X 

X 

MAJ 

X* 

X 

NWA 

SHI 

X 

X 

X 

SCP 

X 

X 

SNZ 

X 

TAT 

X 

X 

X 

ZOB 

X 

Event  1:  13  Sep  80;  OT  21:43:15.7;  4.1S127.5E; 
W  =  29  km;  A/  =  5.5. 

Event  2:  6  May  81;  OT  21:36:06,8;  1.9S  8Q.9W, 

W  =  33  km;  Af  =  6.0. 

Event  3:  4  Mar  81;  OT  21:58:05.9;  38.2N  23.3E; 

H  =  29  km;  Af  =  6.0. 

Event  4:  2  May  81;  OT  16:04:55.6;  36.4N  71, IE; 

H  =  229  km;  M  =  6.3. 

*  Station  MAJ  could  not  be  used  in  slowness  location 
because  estimated  velocity  was  too  hiRh  (36  km/sec).  The 
results  of  these  location  experiments  are  detailed  in 
Table  3. 


TABLE  3 

Event  Locations  on  the  Basis  of  Slowness  and  Azimuth  Observations 


Event 

Number  of 
Station* 

Range 
(degrees  1 

NEir 

Location 

Slownetv*- 

Location 

Azimuth 

Location 

RCC 

SLOW 

(degrees) 

Standard  Error 

Azimuth 

Latitude 

Lonintude 

Latitude 

Longitude 

Latitude 

Longitude 

Latitude 

Longitude 

1 

5 

24-56 

4, IS 

127.5E 

5  1S 

125.9E 

5  7S 

127.7F. 

2.6 

4.6 

3.5 

2 

6 

19-92 

1.9S 

80.9W 

O.OS 

82.  IW 

1.2.S 

79.1  W 

2.0 

17 

1.2 

3 

7(81 

7-S4 

38,2N 

233E 

4  TON 

22  2E 

39.3N 

24.2E 

2.2 

11 

3.6 

4 

8 

30-91 

36  4N 

71  IE 

3.5  ON 

69.3  E 

.34  4N 

69.3E 

3.2 

'it 

8.h 

Tablf  2  data  are  used.  For  the  slownes.s  location,  error  bounds  are  expressed  in  terms  of  the  radius  of  the  65  per  rent 
confidenot  circle  (RCC). 
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TABLE  4 

Fennoscandian  Stations  Used  for  Networking 


Station 

Event  1 

Banda 

Sea 

Event  2 
Ecuador 

Event  3 
Greece 

Event  4 
Afghanistan 

APP 

W 

w 

w 

BER 

W' 

w 

COP 

w 

w 

w 

HFS 

A 

w 

w 

V\' 

KON 

w 

w 

NB2 

W' 

w 

TBV 

w 

W' 

UPP 

w 

w 

w 

SLL 

w 

w 

w 

KEY 

A 

E 

E 

E 

KJF 

A 

E 

E 

E 

KIR 

E 

E 

E 

NUR 

E 

E 

E 

UME 

E 

E 

E 

SOD 

A 

E 

E 

E 

SUF 

E 

E 

E 

ISC  bulletin  reportings  are  used  for  events  listed  in 
Table  2.  E  and  W'  denote  stations  used  in  defining  eastern 
and  western  networks,  respectively.  For  the  Banda  Sea 
event,  a  decomposition  of  Fennoscandian  stations  into 
two  networks  was  not  feasible  due  to  the  limited  number 
of  reporting  stations.  Those  used  for  the  Banda  Sea  event 
are  denoted  with,  A.  Location  results  are  detailed  in 
Table  5. 


TABLE  5 

Summary  of  Networking  Locations 


Event 

Range 

(degrees) 

NEIC 

Location 

Western 

Stations 

Eastern 

Stations 

All  Stations 

Latitude 

Longitude 

Latitude 

Longitude 

Latitude 

Longitude 

Latitude 

Longitude 

1 

94.4 

4.1S 

127.5E 

1.5N 

125.4E 

2 

93.9 

1.9S 

80.9W 

9.8N 

80.7W 

1.39 

81.7W 

O.IN 

80.4W 

3 

23.9 

38.2N 

23.3E 

38.6N 

25.5E 

39.  IN 

22.6E 

39.0N 

23.8E 

4 

42.5 

36.4N 

71. IE 

71.9E 

36.5N 

71.2E 

35.5N 

71.7E 

Event  locations  on  basis  of  fitted  waves  front/slowness  estimation  (Ruud,  1990),  using  the  Fenno- 
■candian  seismograph  network  for  events  and  stations  listed  in  Tables  2  and  4. 


Conclusions 

We  have  here  demonstrated  techniques  for  fast  and  reliable  epicenter  determi¬ 
nation  using  slowness  vector  observations  from  station  networks  and  individual 
three-component  stations.  The  motivation  for  undertaking  the  work  reported  here 
is  the  need  in  many  contexts  of  seismological  data  center  operations  for  fast  and 
reliable  epicenter  locations.  As  demonstrated,  this  can  be  achieved  using  generally 
accessible  seismological  observations,  namely,  parameterized  and/or  original  wave¬ 
form  data.  For  actually  reducing  the  workload  at  various  kinds  of  seismological 
centers,  the  above  location  techniques  must  be  operated  in  an  automated  mode.  In 
the  case  of  networking,  this  entails  that  suites  of  incoming  P  arrival  times  from 
oan\  stations  have  to  be  grouped  in  time  (event  separation)  and  space  (identifica- 
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tiori  of  proper  networks)  before  the  respective  epicenter  locations  can  be  estimated. 
Likewise,  in  the  case  of  waveform  data,  and  given  detection  times,  the  slowness 
vectors  can  easily  be  extracted  together  with  standard  parameters  like  onset  time, 
etc.,  without  analyst  interference.  Now,  given  a  preliminary  location,  arrival  times 
( from  stations  providing  S  observ'ations)  should  be  introduced  for  further  constrain¬ 
ing  the  epicenter  location  (see,  e.g.,  Julian,  1973;  Cassidy,  1989).  Then  it  would  be 
trivial  to  check  whether  P  arrival  reportings  from  other,  conventional  stations  can 
be  associated  with  the  event  in  question  to,  say,  within  10  to  20  sec  of  the  expected 
arrival  time.  This  time  span  is  significantly  shorter  than  it  would  be  were  a 
preliminary  event  location  not  available.  The  P  arrival  window  is  about  20  min, 
including  core  phases. 
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Appendix 

For  estimating  standard  errors  for  azimuth  minimization  we  need  the  0, 
X-derivatives  of  the  G-function.  Consider  a  spherical  triangle  0{‘P,  X),  X) 

and  P(0,  0)  =  North  pole,  where  •P  =  90  —<p  and  <^  =  90  —  0  denotes  colatitude. 
The  following  relation  holds: 

sinla)  =  sin0^)sin(X  —  X)/sinA) 

where  a  is  azimuth,  A  is  epicentral  distance  and  (X  -  X)  is  the  P  angle.  The  relation 
between  A  and  P,  X  is  given  by  the  classic  distance  formula  (Bullen,  1963), 

cos(A)  =  AA  +  B8  +  CC 

where  A  =  sin(S®)cos(X),  B  =  sin(0’)sin(X),  C  =  cos(V^).  Derivatives  with  respect  to 
the  epicenter  parameters  are  obtained  from: 


dG  =  da 


1 

sin(A)cos(a) 


[sinlV?  —  <P)cos(‘P)  dP 


+  sin(<iP)cos(X  —  X)  dX  —  sin(«)cos(A)  tJA] 


and 


-1 

dA  =  — -  [{A  cos(0’)cos(X) 

sm(A) 

+  B  cos(<^)sin(X)  —  sin(<^)cos(V’)]  dP  -  [AS  +  AA]  dX] 

Combining  these  expressions  leads  to  the  derivatives  dGjd^  and  dG/dX.  Finally  the 
axes  of  the  0,  X  confidence  ellipse  can  be  estimated  using  standard  procedures. 
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Exploring  the  upper  crystalline  crust; 

A  joint  interpretation  of  3D  imaging  and 
reflection  profiling  at  the  NORESS  array 

B.O.  Ruud’  and  E.S.  Husebyc’'^ 

‘Dept,  of  geology,  Univ.  of  Oslo,  Norway 
^  Leave  of  absence  from  NORSAR 


Abstract 

P-wavefield  complexities  arc  Ccisily  seen  in  the  NORESS  array  records  de¬ 
spite  its  small  aperture  of  3  km.  The  siting  area  is  well  mapped  geologically  and 
is  partly  ‘sampled’  by  a  2.5  km  seismic  reflection  profile.  Using  this  information 
in  combination  with  tomographic  mapping  (ACH)  on  the  basis  of  NORESS 
Pn-records,  we  have  undertaken  a  structural  mapping  of  the  siting  area  down 
to  ca  3  km.  The  major  result  is  that  the  area  is  relatively  homogeneous  (rhy¬ 
olite/granite)  except  for  the  north-east  quadrant  that  is  dominated  by  an  in¬ 
trusive  gabbro  body.  The  gabbro  intrusion  appears  to  be  roughy  mushroom 
shaped  thinning  towards  the  array  center.  Only  the  upper  profiling  reflector 
(ca  1.1  km)  is  accounted  for  by  this  model  while  deeper  ones  at  ca  1.7  and 
2.5  km  do  not  have  a  counterpart  in  the  3D  results.  The  small  array  aperture 
exclude  mapping  of  more  deep  seated  structural  anomalis. 


Introduction 

Exploration  geophysicists  have  a  preferance  for  seismic  reflection  profiling  a.s  the 
principal  tool  for  remotely  mapping  subsurface  structures  which  otherwise  are  ac¬ 
cessible  only  through  costly  drilling  operations.  On  the  other  hand,  seismologists 
have  a  long  tradition  in  using  wide-angle  reflection  and  refraction  profiling  tech¬ 
niques  in  their  endeavor  to  map  two-dimensionally  the  crust  and  lithosphere  in  ever 
increasing  detail  With  the  advent  of  seismic  tomography  (e.g.,  see  Aki  et  al,  1977) 
a  new  kind  of  mapping  techniques  has  become  very  popular  whenever  observation- 
ally  feasible.  The  reason  for  this  is  naturally  the  possibility  of  obtaining  3D  velocity 
images  of  subsurface  structures. 

Mapping  of  structural  heterogeneities  pose  a  number  of  intrigueing  problems  for 
the  following  reaisons;  petrological  classifications  are  based  mainly  on  mineralogy  and 
textural  criteria,  seismic  reflection  surveys  map  vertical  impedance  contrasts  (first 
order  changes  in  physical  properties),  while  3D  imaging  is  aimed  at  delineating  ve- 
lority  variations  averaged  over  large  volumes.  In  this  context,  the  small  aperture 
array  (3  km)  NORFISS  in  S.-E.  Norway  is  unique;  the  area  is  well  mapped  geologi¬ 
cally,  a  short  (2.5  km)  seismic  reflection  profile  crosses  the  array  (see  Fig.  1)  and  its 
25  vertical  sensors  provide  data  adequate  for  inversion  (imaging)  experiments. 
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Tlie  topic  of  this  paper  is  simply  to  investigate  how  the  reflection  profile  and  3D 
imaging  results  compare  and  to  make  a  joint  interpretation  of  the  two  data  sets. 


Geology  of  the  NORESS  siting  area 

riie  array  is  situated  on  Precambrian  rock  in  southeastern  Norway,  ca  20  kin  east 
of  the  northern  truncation  of  the  Oslo  Graben  structure.  In  Fig.  1  a  small  section 
of  the  geological  maps  (Skalvoll,  1988)  covering  the  NORESS  array  region  is  re¬ 
produced.  The  area  is  dominated  by  granites  and  rhyolites  with  dominant  strike 
southeast-northwest.  Intruded  in  the  granite  and  rhyolites  are  gabbro  bodies  in 
which  physical  properties  (density  and  seismic  velocity)  sharply  contrasts  with  those 
of  the  surrounding  rocks.  The  intrusions  might  either  be  plug-shaped  occasionally 
including  geometrical  ‘mushrooming’  at  the  top  or  siU-shaped  (possibly  dipping  and 
gently  folded).  Subsurface  geometries  can  only  be  hinted  at,  but  as  the  rhyolites  to 
the  southwest  and  northeast  are  connected  to  the  southeast  (outside  the  map  area), 
a  reasonable  hypothesis  is  the  rocks  beneath  the  NORESS  array  have  been  folded 
either  into  a  synform  or  an  antiform. 


The  reflection  profiling  survey 

The  profile  was  shot  in  1983  as  part  of  site  localization  for  the  NORESS  array.  In¬ 
strumentation  used  was  the  Texas  Instrument  DFS  V,  24  geophones  spaced  20  m 
apart  (least  offset  120  m).  Profile  length  is  2  5  km  and  signal  were  generated  with 
dynamite  charges  of  1.4  kg  in  8  m  deep  boreholes  at  40  m  intervals.  During  pro¬ 
cessing  (CMP  stacking)  terrain  corrections  and  crooked  profile  line  corrections  were 
introduced.  For  details  see  Mykkeltveit  (1987). 

The  main  outcome  of  the  profiling  experiment  is  according  to  Mykkeltveit  (see 
also  Fig.  2),  the  existance  of  3  clear  reflectors  within  the  upper  1  sec  TWT  interval, 
namely  (t)  dipping  reflector  at  0.37  sec  TWT  (ca  1.1  km)  at  the  northern  end  of 
the  profile  continuing  up  to  near  surface  at  the  southern  end,  {ii)  dipping  reflector 
at  0.55  sec  TWT  (ca  1.7  km)  in  the  central  part  of  the  profile,  and  (tti)  apparently 
horizontal  reflector  at  0.80  sec  TWT  (ca  2.5  km).  Mykkeltveit  (1987)  tentatively 
interprets  the  uppermost  reflector  as  a  granite/gabbro  interface  dipping  northwards 
at  about  20“’. 


3D  inversion  method  and  results 

The  data  used  for  this  experiment  were  extracted  by  measuring  the  exact  arrival 
time  of  special  features  in  the  seismograms  siwh  as  the  maximum/minimum  ampli¬ 
tude  of  a  peak/ trough  near  the  onset  of  the  signal.  Only  events  with  very  good  SNR 
were  selected  and  the  times  were  determined  with  precision  well  below  the  sampling 
interval  using  natural  cubic  spline  interpolation  between  samples.  By  comparing 
data  estimates  from  several  pair  of  events  at  the  same  locations,  very  high  consis¬ 
tency  was  found,  indicating  a  precision  better  than  2  msec  near  the  onset  of  strong 
signals.  A  total  of  582  relative  arrival  times  from  24  events  with  a  reasonably  good 
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a/iniut  h  (ov.'iaf^o  were  useti  iii  tlio  inversion.  Interestingly,  with  this  kind  of  timing 
precision  array  azimuth  and  apparent  velocities  can  be  estimated  very  accurately, 
that  i',  with  1  2”  and  0.1  0.2  km/sec  respectively,  excluding  bias  effects  from  velocity 
iliniid,  IflS.S;  s('(>  also  Harris,  1!).'<2). 

<il  the  !' t.' ic  assiimirtions  of  tlie  ACil  method,  nec(>.ssary  to  justify  the  use 
■  k,  ;.  ;  ,  -Ji.u-'uiiT.  is  that  tlie  size  (spatial  extension)  of  the  anomalous  struc- 

o.  s  ,  eater  that  the  dominating  wavelength  in  the  data  (c.g.,  for  ACH 

s.n  ‘  .\(ti;SAR  data  wavelengths  are  about  3  km  while  spatial  correlation 
tiic.  c.,  pi.-  kin).  The  dominating  wavelength  in  the  NORESS  records  used  in 
ii.e  '  .  ,  .I'  ll  .;  1  km,  while  the  spatial  correlation  length  of  the  data  residuals 

.  -I.  I  km  I  I;  -  lie-  high  frequency  approximation  inherent  in  travel  time  inver- 
i-'ii  m  .11  icily  inc  lalid  lor  this  ex[)erimenl.  However,  using  principles  from  Fresnel 
diiiiaction  tlmory  wc  may  deline  the ‘width  of  a  ray’  (radius  p)  at  a  distance  D  from 
the  r'-ccivi-r  by  r'Kiuiring  that  the  difference  in  distance  to  the  diffraction  point  in 
ilic  aavefront  to  be  hess  than  a  quarter  of  the  wavelength  A  so  that  they  interfere 
coiotructis'ely  (Nolet,  11)87,  p.  10).  For  a  source  at  infinite  distance  we  then  have 

7^2  +  £)•;_  D  =  A/.f  (1) 

.8ol\  ing  for  p  wo  get 

P=^{WTd^  (2) 

l  or  in'innco.  for  A  —  I)  —  1  km  we  get  p  =  0.75  km  which  gives  a  clue  to  the 
'r'  .-,oIving  power’  of  tiie  data  at  a  depth  of  1  km.  This  reasoning  was  used  when 
-I'-cidiiig  the  size  of  the  blocks  used  in  the  tomographic  experiment. 

.Aunt  h('r  basic  assumption  of  the  method  is  that  the  anomalies  should  be  confined 
1 1 1  the  interior  of  the  model  so  that  the  wave  is  plane  as  it  enters  into  the  model  at  the 
bottom  interface.  This  assumption  is  in  most  cases  not  strictly  valid,  but  in  practice 
the  model  is  made  a  little  larger  than  the  volume  we  expect  to  be  well  resolved  by 
the  data.  Any  velocity  anomalies  outside  but  dose  to  the  model  will  then  appear 
piim  lily  in  the  bottom  layer.  Larger  scale  (larger  that  the  model)  anomalies  far 
from  the  array  could  refract  the  wave  such  that  it  arrive  from  a  direction  somewhat 
lifferent  from  the  one  expected  from  bulletin  information,  but  still  mainly  as  a  plane 
wave.  To  account  for  this  we  define  travel  time  anomalies  relative  to  the  best  fitting 
plane  wave  and  not  relative  to  expected  azimuth  and  slowness.  Denoting  arrival 
time  of  the  i-th  event  at  the  j-th  sensor  as  LjP  the  average  arrival  time  as  t,-,  the 
slowness  vector  as  u,,  and  the  position  of  the  j-th  sensor  relative  to  the  array  centre 
as  r_,,  the  time  residuals  are  defined  as 

=  l,j~  (i,+  Tj-  u,)  (3) 

le-re  H  and  u,  are  lea.st  square  <-stimates  for  each  event.  Elevation  corrections  were 
'■  l,,r  (although  .lot  included  in  the  equation  above). 

In  <  -!  Hast  to  the  original  ACH  method  all  event  parameters  were  retained  as 
unknown  through  the  inversion  proec  s.  The  estimated  event  slownesses  are  used 
1)1  r.iy-tracing  in  the  initial  lateral  homogeneous  model  in  order  to  compute  the 
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It'ngUi  of  tlu'  ray  soginciits  diji;  (cvonl  i,  sensor  j,  block  k).  'I'lic  system  of  equations 
to  invert  is  tlien 

Af.j  =  At,  +  Tj  ■  All,  +  Asa  (4) 

k 

where  At,  and  Au,  are  correction  terms  to  tlie  event  parameters  and  As  is  the 
slowness  anomalies.  'J'his  e.xtcnded  system  was  solved  using  the  metliod  of  parameter 
separation  (Spencer  and  Cnbbins,  1980)  which  reduce  the  size  of  the  system  to 
the  same  ;us  for  the  original  AClf  method.  Stochastic  inversion  was  used  to  dam[) 

*  unre.solved  model  parameters. 

The  obtained  slowness  anomalies  are  presented  in  Fig.  3  together  with  standard 
deviations.  The  variations  in  standard  deviation  depends  mainly  on  model  resolution 

^  which  is  very  good  in  the  central  blocks  in  each  layer  (about  0.9),  but  decre.asing 

towards  the  peripheral  blocks.  As  usual  with  ACll  typc  of  inversion  horizontal 
velocity  layering  can  not  be  resolved  because  relative  time  residuals  are  used.  The 
RMS  value  of  the  residuals  were  5.5  msec  which  is  reduced  to  3.7  msec  by  the 
inversion.  In  terms  of  variances,  the  model  explains  53%  of  the  observed  residuals. 
The  lateral  dimensions  of  the  model  were  increased  with  depth  in  order  to  assure  that 
all  rays  are  within  the  model  volume  when  entering  the  bottom  layer  from  below 
(the  incidence  angle  of  most  rays  is  less  than  47°,  corresponding  to  Fn  velocity, 
ca  8.2  km/sec)  and  the  block  size  were  also  increased  according  to  the  expected 
<l‘'crf'ase  in  resolving  power  with  depth. 

Regarding  the  estimated  slowness  anomalies,  we  see  that  the  most  significant 
anomalies,  i.o.,  anomalies  greater  than  the  estimation  error,  are  confined  to  the 
uppermost  layer  (depth  range:  0.0-1. 0  km)  and  that  there  is  a  dominace  of  high  ve¬ 
locities  to  the  north-east  in  this  layer.  In  the  second  layer  (depth  range:  1. 0-2.5  km) 
anomalies  are  generally  much  smaller,  and  we  note  that  the  strongest  anomalies  are 
confined  to  the  peripheral  blocks  where  estimation  errors  also  are  large.  For  the 
third  layer  (depth  range:  2. 5-5.0  km)  slowness  anomalies  are  so  weak  that  this  layer 
is  considered  essentially  homogeneous.  Thus  there  seems  to  be  no  reason  to  add 
more  layers  to  the  model  which  besides  would  give  a  depth  range  far  exceeding  the 
array  aperture  of  3  km. 

Discussion 

*  From  the  geology  of  the  NORESS  siting  area  our  working  hypothesis  in  discussing 
the  seismic  results  is  that  the  observed  reflectors/anomalies  are  tied  to  differences  in 
physical  properties  between  granite/ rhyolite  and  gabbro  at  least  for  the  uppermost 

4  one  at  ca  1.1  km  (layer  1).  On  this  account  we  interpret  this  reflector  as  the  underside 

of  an  intruding  gabbro  body.  From  the  3D  results  the  inferred  body  appears  to  be 
steeply  dipping  northeast-  wards,  while  the  profile  gives  a  more  moderate  northward 
dip.  Note  that  these  results  are  not  in  conflict  because  the  profile  is  not  perpendicular 
to  the  strike.  According  to  the  geological  map  a  reflection  from  the  top  of  the  gabbro 
intrusion  should  possibly  be  present  in  the  profile  section.  However,  the  quality  of 
the  upper  2-300  meters  of  the  section  is  poor  due  to  interference  from  the  dominant 
direct  waves. 
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'I'lio  wosteni  and  soutlicni  parts  of  tlic>  N()l{  JCSS  siting  area  appears  to  be  rather 
lioinogcnoous  in  layer  2.  'I'lie  31)  results  imply  that  granite/rhyolite  rocks  dom¬ 
inate  here  in  view  of  the  relative  low  velocities  and  this  hy])othesis  is  consistent 
with  the  local  geology.  On  the  other  hand,  relative  high  velocities  j)ersist  in  the 
north-east,  and  hence  tied  to  a  thickening  <}f  tin-  inferred  gabbro  body.  This  in¬ 
terpretation  impli''s  that  the  1.7  and  2.5  km  rellcctors  arc  unlikely  to  be  associated 
with  a  granite-rliyolite/gabbro  sandwiching  structure.  In  short,  the  origin  of  the  two 
deeper  reflectors  remain  obscure  with  the  data  at  hand  and  we  can  only  regret  that 
there  is  no  profding  perpendicular  to  the  dominant  strike  direction.  However,  we 
note  that  observations  of  seismic  reflectors  in  the  crystalline  crust  arc  widespread 
although  their  geological  counterpart  often  remain  diffuse  (e.g.,  see  Larsson  and 
llusebye,  19S9). 

In  the  bottom  layer  the  3D  anomabes  arc  weak  and  besides  cannot  be  given  any 
geological  significance.  To  summerize  the  discussion  so  far;  the  3D  anomabes  are 
interpreted  in  terms  of  a  gabbro  body  in  the  north-east  quadrant  of  the  array  siting 
area  thickening  towards  the  north-east.  The  profiling  data  seemingly  just  sample 
the  flank  of  a  mushroom-bke  intrusive  body. 

In  the  profibng  results  there  is  indication  of  a  weak  reflector  at  ca  6  sec  TWT 
(Mykkeltveit,  1987  -  not  shown  in  Fig.  2)  or  equivalently  at  a  depth  of  ca  18  km.  This 
depth  is  ca  6  times  larger  than  the  NORESS  aperture  for  which  there  is  no  resolution 
in  the  3D  analysis.  However,  in  3-component  wavefield  studies  (e.g.,  see  Christof- 
fersson  et  ah,  1988)  using  NORESS  data  occasionally  SV-wavelets  are  observed  at 
ca  2  sec  after  the  first  P-onset  and  the  corresponding  ‘conversion’  boundary  would 
be  about  18  km.  Also,  gravity  studies  over  the  Oslo  Graben  (Wessel  and  Husebye, 
1987)  implies  the  existence  of  an  intrusive,  magmatic  body  at  this  depth  range. 

Finally,  to  what  extent  should  we  really  expect  overlapping  results  using  rather 
different  mapping  techniques?  Reflection  profiling  essentiaUy  maps  discontinueties  in 
acoustic  impedance  along  near  horizontal  planes  using  high  frequency  (100-250  Hz) 
reflected  waves  and  hence  good  depth  resolution  (within  50  m).  With  3D  imaging, 
using  direct  waves  from  distant  sources  in  the  frequency  range  2-10  Hz,  we  map 
relative  long  wavelength  (of  the  order  of  1  km)  velocity  variations.  The  inherent 
smoothing  involved  would  naturally  reduce  velocity  anomabes  compared  to  those 
for  web  separated  gabbro/granit  bodies.  In  the  NORESS  siting  area  it  is  not  easy 
to  reconcile  reflection  and  3D  imaging  results  -  and  in  this  respect  being  hampered 
by  very  sparse  reflection  coverage  -  most  needed  is  a  profile  perpendicular  to  the 
NW  strike  direction. 


Conclusions 

In  this  study  we  have  demonstrated  that  even  though  seismic  profiling  and  3D 
inversion  results  seemingly  are  very  different  the  two  methods  actually  provide  com¬ 
plementary  information  to  the  geological  interpretation.  While  the  reflection  data 
maps  rock  boundaries  with  high  re.solution,  the  3D  imaging  data  provide  information 
on  velocity  variations  helping  to  classify  the  structural  volumes  and  to  indicate  their 
spatial  extension.  Further  profibng  remains  desirable  for  a  more  detailed  mapping 
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of  the  structures  and  for  a  better  comparison  of  the  two  methods. 

The  seismic  results  together  with  local  geology  information  strongly  imply  that 
the  north-eastern  part  of  the  NORESS  siting  area  is  underlain  by  a  shallow  gabbro 
intrusion  possibly  mushroom  shaped  and  thickening  towards  the  north-east.  We 
intuitively  expect  that  this  body  would  influence  the  wavefield  as  observed  across 
the  array  and  this  appears  indeed  to  be  the  case.  For  example,  Christofiersson 
(pers.  communication)  reports  a  buckling  of  the  P  wavefront  as  reconstructed  inde¬ 
pendently  from  estimated  slowness  vectors  using  3-component  stations  within  the 
array.  Ultimately,  we  plan  to  undertake  3D  synthetic  seismogram  analysi  to  check 
to  what  extent  the  gabbro  body  beneath  NORESS  can  account  for  observed  P-wave 
train  complexities  as  observed  across  the  array. 
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Figure  captions 


Figure  1;  Geological  map  over  the  NORESS  array  and  its  vicinity  (array  center  at 
60.74°N  and  il.54°E).  The  array  seisinonicters  used  in  the  inversion  experiment  are 
marked  with  black  dots  and  the  reflection  profile  line  with  a  ‘P’.  Typical  density  and 
P-velocity  values  for  granite  and  gabbro  are  (2.66  g/cm^,  5.9  kin/s)  and  (2.99  g/cm^, 
6.7  km/s)  respectively.  The  physical  properties  of  rhyolite  are  expected  to  be  close 
to  those  of  granite. 

Figure  2:  Section  of  tlie  reflection  profile  showing  the  upper  1  sec  TWT  (down  to 
ca  3  km).  TWT  is  two-way-traveltime  in  seconds.  The  three  dominant  reflectors 
referred  to  in  the  text  are  marked  Ri,  R2  and  R3.  Tlie  traces  have  been  filtered  in 
the  100-250  Hz  passband. 

Figure  3:  Slowness  anomalies  obtained  from  inversion.  The  lower  numbers  in 
each  block  are  standard  deviations.  Both  anomalies  and  uncertainties  are  given 
in  msec/km  relative  to  the  reference  layer  velocities  (5.8,  5.9  and  6.0  km/s).  A 
slowness  anomalie  of  -t-3  msec/km  would  thus  correspond  to  a  velocity  change  of  ca 
-0.1  km/s.  Labels  H  and  L  indicates  high  and  low  velocity  regions.  Only  significant 
anomalies,  i.e.,  greater  than  standard  deviation  are  hatched.  Note  the  differences 
in  length  scale  for  each  layer  (indicated  in  the  figure).  The  lateral  block  sizes  are 
0.6  km  for  layer  1,  0.9  km  for  layer  2  and  1.5  km  for  layer  3. 


25 


1 

o 

1 

C\J 

DC^. 

1 

CD 

00 

q 

d 

o 

d 

d 

lAAi 

27 

0.7  0.\  -0.5  0.6  0.3 

2.4  \(.7/  \.\\  l.d  1.0  1.1  1.3 


-1.0  O.lv  Xo.9/  /l.Aj  0.1  0.1  -0.2 


.OV  1.0  1.1  1.4 


2.5  1.3Y  1.] 

V  y 


-0.2  0 . 3  A  .  ^  I  1 .  tN  -0 . 7  -0.4  -0.7 

2.5  1.3  1.0  1.2  3.1 

0.0  \  .X<r\  0.4  0.6  0.6 

l/^xV-2.>l.l  1.1  1.3  3.0 


Figurf*  3:  Slowness  anomalies  obtained  from 
inversion.  The  lower  numbers  in  each  block 
are  standard  deviations.  Both  anomabes  and 
uncertainties  are  given  in  msec/km  relative 
to  the  reference  layer  velocities  (5.8,  5.9  and 
6.0  km/s).  Labels  H  and  L  indicates  high 
and  low  velocity  regions.  Only  significant 
anomalies,  i.e.,  greater  than  standard  devi¬ 
ation  are  hatched.  Note  the  differences  in 
length  scale  for  each  layer  (indicated  in  the 
figure).  The  lateral  block  sizes  are  0.6  km 
for  layer  1,  0.9  km  for  layer  2  and  1.5  km  for 
layer  3. 
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ABSTRACT 

Pioneering  work  on  mapping  the  Scandinavian  crust  commenced  in  the  early  1960s  and  since 
then  numerous  profiling  surveys  have  been  undertaken,  particularly  as  part  of  the  on-going 
EUGENO-S  project.  However,  the  most  significant  contribution  to  mapping  crustal  structural 
details  came  from  the  MA^  Mobil  Search  cruises  in  the  Skagerrak  Sea  and  off  the  West 
coast  of  Norway  (reflection  profiling  16  sec  TWT).  In  short,  all  past  and  present  crustal 
profiling  results  have  been  integrated  to  produce  detailed  maps  of  Moho  depths  and  crustal 
thicknesses  for  the  S.  Scandinavia  The  thinnest  crust  is  found  in  the  North  Sea  and 
Skagerrak  (approximatly  20  km),  while  east-central  Sweden  feature  very  thick  crust  (app.  50 
km).  Other  interesting  features  are  the  apparent  correlation  between  crustal  thiiming  and 
sedimentation  /  subsidence,  documented  magmatic  activity,  earthquake  occurrences  and  the 
tectonic  age  of  the  crust.  Moho  depth  and  the  crustal  thicknesses  clearly  mirror  the  tectonic 
evolution  and  the  present  structural  features  of  the  region  investigated. 
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1. 


INTRODUCTION 


Crustal  studies  became  popular  with  seismologists  some  3  decades  ago.  In  the  early 
sixties  the  prime  motivations  were  to  map  Moho  depths,  deduce  whether  the  number  of 

▲ 

crustal  layers  should  be  2,  3  or  4,  and  investigate  to  what  extent  velocity  reversals  took 
place.  Now  as  before  the  most  popular  and  also  the  most  precisattemte  tool  for  such  studies 
are  tied  to  seismic  surveys,  i.e.  refraction  and/or  wide  angle  reflection  profiling.  However, 
during  the  last  decade  improved  instrumentation  in  combination  with  sophisticated  forward 
and  inverse  modelling  has  permitted  more  detailed  mapping  of  crustal  structures.  Though  1-D 
line  profiling  is  limited  to  identifying  gently  undulating  velocity  contrasts  in  the  crust  and 
below  Moho,  and  naturally  depths  to  majore  structural  boundaries.  Significant  improvements 
in  mapping  came  with  the  advent  of  seismic  tomography  (e.g.,  see  Aki  et  al.,  1977)  and  the 
adaption  to  crustal  studies  of  the  seismic  reflection  techniques  used  in  oil  prospecting,  (e.g. 
the  USA  COCORP,  and  the  UK  BIRP  programs).  The  use  of  these  mapping  techniques 
provide  vastly  improved  structural  images  of  crustal  heterogeneties,  but  their  inteipretation  in 
a  geological  context  remains  non-unique.  For  example,  the  abundance  of  crustal  reflectors  or 
laminae  observed  in  deep  crustal  reflection  surveys  can  be  explained  in  terms  of  petrological, 
geochemical,  hydrostatic  or  magmatic  discontinuities  or  simply  low-angle  faults  (Larsson  and 
Husebye,  this  volume). 

The  crustal  study  developments  in  Scandinavia  have  a  clear  parallel  to  those  taking 
place  elsewhere;  it  started  with  refraction  profiling  in  the  early  sixties,  extended  to  wide-  * 

angle  profiling  and  tomography,  and  recently  culminated  with  approximatly  10.000  km  of 
deep  seismic  reflection  data  (16  sec  TWT)  collected  during  M/V  Mobil  Search  cruises  in 
Norwegian  off-shore  areas  in  1987.  Our  part  of  this  project  was  the  Skagerrak  Sea  survey, 
totaling  1730  km  of  intersecting  profiling  lines.  The  shotpoint  map  is  shown  in  Fig.  1 
(Husebye  et  al.,  1988).  The  crustal  information  extracted  from  these  marine  data  are  the 
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main  topic  of  this  presentation.  Fortunately,  it  is  fairly  easy  to  integrate  the  Skagerrak  survey 
into  seismic  profiling  eftorts  of  the  past,  thus  the  study  was  extended  to  comprise  the  whole 
of  south  Scandinavia.  As  we  consider  it  important  to  attempt  to  interpret  geophysical  results 
in  a  geological  context,  we  therefore  start  with  a  brief  overview  of  the  geological  evolution 
of  not  only  south  Scandinavia  but  the  whole  of  Fennoscandia. 


2.  TECTONIC  EVOLUTION  OF  FENNOSCANDIA  WITH  EMPHASIS  ON  THE 
SOUTH  SCANDINAVIAN  AREA 


A  structural  map  of  south  Scandinavia  is  shown  in  Fig.  2.  Major  stmctural  features 
and  evolutionary  details  are  as  follow; 

2.1  BALTIC  SHIELD  STRUCTURES 

The  Baltic  Shield,  comprising  mainland  Norway,  Sweden,  Finland,  the  Kola  peninsula 
and  USSR  Karelia,  was  mainly  formed  or  accreted  during  four  orogens  (enhanced  in  Table 
I)  essentially  confined  to  the  time  period  3100  to  1500  Ma.  Later  orogens  have  largely 
reworked  the  existing  crust,  adding  little  if  any  new  material  to  the  craton.  The  shield  is 
roughly  zoned,  in  the  sense  that  the  oldest  cnist  is  found  in  NE  Finland  and  the  Kola 
peninsula  and  the  youngest  in  southern  Norway.  ^Tiile  the  oldest  crust  stems  from  the 
Saamian  orogen  (3100  -  2900  Ma),  most  of  the  Baltic  Shield  was  formed  during  the 
following  three  orogens;  the  Lopian  (2900  -  2600  Ma),  the  Svecofennian  (2000  -  1750  Ma) 
and  the  Gothian  (1750  -  15(X)  Ma). 

The  Lopian  and  Svecofennian  orogens  had  a  northwest-southeasterly  front  while  the 
Gothian  was  more  north-south  (Berthelsen  and  Marker,  1986).  According  to  Poorter  (1981) 
and  Steam  and  Piper  (1984),  the  Baltic  Shield  was  rotated  clockwise  with  respect  to 
Laurentia  during  the  interval  1190  to  1050  Ma.  This  also  account  for  the  complexity  of  the 
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Gothian  orogen  (Gadl  and  Gorbatschev,  1987). 

The  Protogine  zone  (PZ  in  Fig.  2)  in  south-central  Sweden  was  formed  by  episodic 
faulting  accompanied  by  magmatic  intrusions.  Main  phases  of  activity  took  place  between 
1650  -  900  Ma.  Today  it  constitutes  a  sharp  boundary  between  the  highly  metamorphic  and 
foliated  SW  Scandinavia  (reworked  during  the  Greenvillian  orogen)  and  the  well  preserved 
Svecofennian  terrains  to  the  east.  The  effect  of  the  later  Sveconorwegian  orogen  was  a 
thorough  reworking  of  the  crust  as  far  east  as  the  Protogine  zone. 

During  the  Caledonian  orogen  (600  -  400  Ma)  large  areas  of  the  precambrian  crust  of 
west  Scandinavia  became  concealed  beneath  nappes,  consisting  of  disconnected  precambrian 
basement  fragments  with  younger  sediments  ontop.  This  orogen  affected  the  precambrian 
basement  progressively  in  a  westward  direction  by  reworking  and  foliation.  The  Oslo  Graben 
(OG)  is  an  intracontinental  rift,  cutting  through  the  younger  parts  of  the  Baltic  Shield  and  is 
dominated  by  large  amounts  of  volcanic  material,  dating  from  305  Ma  in  the  south  to  245 
Ma  in  the  north  (Neumann  et  al.,  1986). 

We  have  only  given  a  rough  outline  of  the  geological  evolution  of  the  Baltic  Shield 
for  the  simple  reason  that  such  information  is  difficult  to  correlate  with  seismic/geophysical 
data.  The  fundamental  problem  is  that  past  deformations  are  far  better  preserved  geologically 
and  geochemically  than  in  geophysical  anomalies  refreciing  the  pieseui  tectonic  environment. 

2.2  DENMARK  AND  ADJACENT  SEAS. 

Denmark  and  adjacent  seas  belong  to  the  mobile  part  of  paleo-Europe,  and  as  such 
appear  to  be  younger  than  the  stable  shield  areas  to  the  north.  How  and  when  these  areas 
became  attached  to  each  other,  remains  somewhat  obscure.  Berthelsen  (1987)  has  suggested 
that  a  continental-continental  collision  took  place.  Denmark  and  adjacent  seas  have  been 
subject  to  more  recent  and  severe  deformations  than  the  shield  areas  to  the  north.  This  is 
demonstrated  in  a  Skagerrak  basin  subsidence  analysis  by  Pedersen  et  al.  (this  volume),  who 
also  reconstructed  crustal  thicknesses  to  precambrian  times.  From  Fig.  2  we  see  that  this  area 
is  dominated  by  three  main  features,  namely  the  Fennoscandian  Border  Zone  (FBZ),  the 
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Norwegian-Danish  Basin  (NDB),  and  the  Ringk0bing  Fyn  High  (RFH).  The  FBZ  is  a  NW  to 
NNW  tiending  fault  zone  separating  the  Shield  and  the  European  platform.  Even  though  the 
FBZ  has  been  repeatedly  reactivated  since  precambrian  times,  the  movements  along  the  zone 
appear  modest  (Liboriussen  et  al.,  1987;  Komfaelt  and  Larsson,  1987).  The  sediment 
thicknesses  in  the  NDB  are  large,  exceeding  10  km  in  NW  Jutland.  Subsidence  probably 
started  in  early  Paleozoicum  but  the  main  pulse  took  place  during  the  Triassic.  The 
sediments  here  are  dominated  by  thick  Zechstein  evaporites.  Subsequent  halokinetic 
movements  resulted  in  domes  and  diapirs  in  the  NW  part  of  the  basin.  The  RFH  is  an  area 
of  elevated  precambrian  basement  and  is  subdivided  into  several  blocks  by  N-S  trending 
grabens  such  as  the  Horn  Graben.  RFH  appears  to  have  been  elevated  in  the  Carboniferous 
while  taphrogenesis  has  been  related  to  the  late  Carboniferous  -  early  Permian  Hercynian 
movements  (Ramberg  and  Spjeldnaes,  1978;  EUGENO-S  Working  Group,  1988).  A  brief 
overview  of  the  evolution  of  this  area  is  given  in  Table  II.  A  very  extensive  literature  exists 
on  the  evolution  of  Denmark  and  adjacent  seas,  and  key  references  here  are  Fluh  and 
Berthelsen,  1986;  Michelsen,  1976;  Ziegler,  1981  and  1982  and  Pegrum  1984a  and  1984b. 


3.  SEISMIC  PROFILING  DATA  AND  THEIR  ANALYSIS 


Initially,  the  intention  was  to  restrict  the  crustal  analysis  to  the  M/V  Mobil  Search 
profiling  lines  in  the  Skagerrak  Sea  (Fig.  1)  where  reflection  records  are  of  exceptionally 
high  quality.  The  opportunity  for  tying-in  the  Skagerrak  crustal  results  with  those  obtained 
from  nearby  land  profiling  surveys  made  it  natural  to  extend  the  scope  of  our  crustal 
synthesis  to  the  whole  of  S.  Scandinavia.  This  offer  an  unique  opportunity  to  examine  crustal 
variations  across  vastly  different  tectonic  provinces.  The  Skagerrak  survey  is  dealt  with  in 
detail  below  while  for  the  other  profiling  surveys  (shown  in  Fig.  3)  we  have  used  the 
relevant  results  as  summarized  in  reports  and  publications  (see  Table  IV). 
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3.1  THE  MA^  MOBILE  SEARCH  SURVEY  -  DATA  AND  ANALYSIS 

During  the  Mobile  Search  caiise  of  one  month  in  Winter  1987  (Husebye  and  Ro, 
1987;  Husebye  et  al.,  1988)  1730  km  of  seismic  reflection  profiling  was  accomplished.  A 
dual  sampling  rate  of  4  and  8  msec  was  used,  and  record  length  was  16  sec  TWT  giving  an 
effective  penetration  depth  of  50-60  km.  On-board  data  processing  at  8  msec  and  the 
resulting  seismic  sections  constitute  the  basis  for  our  cmstal  analysis  in  the  area.  Data 
quality  is  exceptionally  high,  and  may  be  further  enhanced  by  reprocessing  using  more 
refined  techniques  than  permitted  on-board.  Such  a  project  is  now  being  completed  with 
financial  support  from  oil  companies. 

The  first  step  in  analysis/interpretation  of  the  individual  profiles  was  to  identify  major 
reflectors  in  the  sedimentary  column  and  crystalline  crust.  Guided  by  other  surveys  and 
works  in  the  area,  notably  that  of  Haatvedt  and  Nerby  (1984)  and  bore-hole  log  information 
from  the  Danish  Geological  Survey,  reflectors  were  tied  to  the  following  horizons:  Moho, 
Precambrian  basement,  top  Permian,  base  and  top  Zechstein,  inter-Jurassic  erosional  plane, 
base  and  top  Chalk  and  base  Tertiary  and  Quartemary.  Depth  conversion  was  achieved  using 
standard  techniques  and  velocities  as  listed  in  Table  III.  The  results  of  the  analysis  are 
illustrated  in  Fig.  4  and  5  for  OG-2  and  12  and  5,7  and  13  respectively.  On  the  basis  of 
individual  profiling  results,  maps  of  Moho  depths  and  thickness  of  crystalline  crust  have 
been  constructed  and  are  displayed  in  Figs.  6  and  7. 

A  few  clarifying  comments  to  these  maps  apply;  the  quality  and  the  reflector 
sharpness  varies  from  line  to  line  and  along  individual  lines.  In  the  southernmost  profiles 
halokinetic  make  it  difficult  to  extract  information  from  reflectors  below  Zechstein;  this 
problem  is  most  pronounced  on  OG-4.  In  north  Skagerrak  a  combination  of  thin  sediments 
and  a  hard  bottom  with  a  rough  relief  make  discrimination  between  the  paleozoic  sediments 
and  basement  often  non-unique.  In  other  words,  the  exact  position  of  the  southward 
continuation  of  the  Oslo  Graben  (OG),  that  is  the  Oslo  Rift  (OR),  is  somewhat  diffuse.  On 
some  profiles  Moho  can  be  tied  to  a  sharp  and  well  defined  reflector,  characterized  by  2 
cycles  of  high  amplitude.  On  others,  there  are  bands  or  several  refleaors  in  the  appropriate 
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depth  range.  This  is  not  tied  only  to  line  location  but  also  to  profile  strike  direction.  The 
sharp  reflectors,  like  those  on  OG-13.  have  been  used  as  markers  for  identifying  Moho  in 
reflector  bands  at  intersections  with  other  lines,  such  as  OG-12.  A  basic  problem  here  is 
naturally  that  die  Moho  discontinuity  is  not  uniformly  sharp  over  the  entire  area.  However, 
the  correlation  between  OG-12  and  OG-13  indicates  that  Moho  should  not  be  defined  as  the 
bottom  of  a  general  and  broad  reflective  zone  which  seems  to  be  the  tendency  in  many 
recent  works  in  deep  seismic  reflection  profiling  (e.g.,  see  Matthews  and  Cheadle,  1986). 

To  summarize,  data  processing  tied  to  8  msec  sampling  is  a  bit  coarse  for  the 
sedimentary  column  but  considered  adequate  for  identifying  majore  discontinuities  and 
detailing  velocity  distributions.  The  final  outcome,  that  is  Moho  depth  estimates,  would  as 
mentioned  depend  on  choice  of  reflector(s)  taken  to  be  representative  for  this  prime 
discontinuity.  It  is  not  always  sharp  (Barton,  1986)  but  by  harmonizing  reflector  "pickings" 
at  line  intersections  a  reasonable  consistent  Moho  isopac  map  has  been  obtained.  Errors  here 
are  considered  unlikely  to  exceed  0.5-1  sec  TWT  or  equivalently  that  Moho  depth  estimates 
are  correct  within  ±3  km,  while  the  relative  error  would  be  less.  In  the  south  Moho 
reflectors  are  more  diffuse  and  we  may  even  have  a  more  bumpy  Moho  across  the 
Fjerritslev  Fault  Zone  (FFZ),  and  the  depth  errors  could  be  slightly  larger  here.  During 
reprocessing  special  attention  would  be  given  to  this  and  similar  problems.  Finally, 
upper/lower  crust  transition  is  not  a  first  order  discontinuity  in  Skagerrak,  as  reflectorfs)  at 
appropriate  depth  range  have  not  been  identified.  However,  from  extensive  laminae  analysis, 
the  thickness  of  the  upper  brittle  and  lower  ductile  parts  of  the  crust  appears  to  be  roughly 
the  same,  that  is  50  ±  10  per  cent  (Larsson  and  Huseby,  this  volume). 


3.2  SEISMIC  REFRACTION  AND  WIDE  ANGLE  PROFILING  SURVEYS 

In  the  past  numerous  seismic  profiling  surveys  have  been  conducted  over  various 
parts  of  south  Scandinavia  -  many  of  those  bordering  on  the  Mobile  Search  profiling  lines  in 
Skagerrak.  Details  on  these  profiles  are  given  in  Fig.  3  and  in  Table  IV.  Note  that  some  of 
the  oldest  profiling  records  (shot  prior  to  1971)  have  been  reinterpreted  by  Kvsema  (1984), 
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these  are  profiles  I,  K  and  L.  The  Moho  depth  estimates  from  all  these  surveys  plus  those 
from  Skagerrak  has  been  integrated,  the  final  result  being  the  south  Scandinavian  Moho 
depth  and  crustal  thickness  maps  in  Figs.  6  and  7.  Noteworthy,  the  Skagerrak  results  blend 
well  with  those  of  bordering  profiles  in  Norway,  Sweden  and  Denmark  -  differences  being 
within  ±2km.  Naturally,  in  drawing  the  Moho  depth  map  some  cosmetic  smoothing  has  been 
introduced  in  part  justefied  by  covering  gaps  between  adjacent  profile  lines.  Our  mapping 
should  be  considered  an  evolutionary  process  i.e.,  it  is  updated  as  new  profiling  results 
become  available. 

4.  RESULTS 

The  Moho  depth  and  the  cmstal  thickness  may  in  comparision  to  the  structural 
features  indicate  the  relative  correlation  between  these  features  and  thus  provide  a  clue  to  the 
evolution  of  the  region  under  investigation. 

The  main  features  of  the  Moho  depth  map  (shown  in  Fig.  6)  is  a  pronounced  Moho 
elevation  in  Skagerrak,  Jutland  and  off-coast  west  Norway,  which  geological  counterparts  are 
sedimentary  basins,  including  the  NDB.  This  is  expected  as  basin  formation  is  a  direct 
consequence  of  extension  processes  or  crustal  thinning  (e.g.,  see  Pedersen  et  al.,  this 
volume).  In  particular,  there  is  a  marked  Moho  elevation  of  approximately  5km  associated 
with  the  OR/OG  system,  being  somewhat  broader  than  the  corresponding  surface  structural 
features.  The  crust  underneath  the  Caledonides  of  west  Norway  is  relatively  thick,  while  the 
Western  Gneiss  Region  is  thin.  There  is  a  similar  feature  of  thicker  crust,  also  related  to  the 
Caledonides  on  the  British  side  ofViking  Graben/Central  Graben  (VG/CG),  but  a  correlation 
is  speculative  (S.  Klemperer  pers.  comm.).  Pronounced  crustal  thinning  take  place  off  S. 
Norway  and  SW.  Sweden.  In  the  laner  case  the  trends  of  the  FFZ/FBZ  seem  to  be 
dominating.  The  RFH  stands  out  sharply  as  a  thick  crustal  block.  The  deepest  Moho  is 
found  in  southeast  Sweden  with  depths  around  50  km.  The  Moho  trends  in  Skagerrak  seem 
to  be  influenced  mainly  by  the  OG/OR  and  the  NDB.  The  FFZ/FBZ  seem  to  be  of  less 


importance  in  this  area,  but  its  more  easterly  segments  (in  the  Kattegat  area)  appear  to  have 
a  counterpart  in  Moho  undulation.  The  oldest  part  of  the  Baltic  Shield,  the  area  east  of 
Protogine  Zone(PZ)  or  the  Sveconorwegian  FrontfSNF),  is  relatively  uniform  with  depths 
exceding  40  km  and  only  minor  gradients. 

The  thickness  of  the  crystalline  crust  is  displayed  in  Fig.  7.  The  thinning  in 
Skagerrak  and  Denmark  is  pronounced  with  a  minimum  in  NDB  of  about  18  km.  Again  the 
pattern  is  complex  with  trend  directions  coinciding/following  FFB/FFZ,  Swedish  west  coast, 
the  OR/OG  and  fmally  one  paralleling  the  SE  coast  of  Norway.  The  RFH  is  even  more 
pronounced  than  in  the  Moho  depth  map,  and  an  off-shore  westward  extending  "tongue"  of 
the  Caledonide  province  is  clear,  fhe  extreme  localized  thinning  near  Bornholm  (22  km) 
coincides  with  a  southward  shift  in  FEZ.  Areas  not  affected  by  post-Cambrian  subsidence 
and  rifting  coincides  generally  with  the  area  where  the  crust  is  thicker  than  36  km. 


5.  OTHER  GEOPHYSICAL  AND  GEOLOGICAL  OBSERVATIONS 

It  is  prudent  to  consider/discuss  our  cmstal  results  in  the  framework  of  other 
geophysical  and  geological  observations,  and  a  brief  summary  of  these  are  given  below. 


5.1  VOLCANIC  ACTIVITY  IN  SOUTH  SCANDINAVIA:  CARBONIFEROUS  TO 

Present 

The  western  part  of  S.  Scandinavia  has  experienced  volcanic  activity  as  shown  in 
Fig.  8.  Paleozoic  magmatism  is  mainly  related  to  the  formation  of  th  Oslo  Rift  (OR)  and  the 
Fjerritslev  Fault  Zone  (FF^.  Manifestions  on  land  are  dykes  in  S.  Scandia  (Bergstr0m, 
1982;  Klingspor,  1976)  and  of  dykes  and  magma  outcrops  in  the  Oslo  Graben  (OG) 
(Neumann  et  al.,  1986;  Ramberg  and  SpjeldntBS,  1978;  Ramberg  and  Larsen,  1978).  In 
Denmark  and  Danish  North  Sea  magma  of  this  age  have  been  located  in  some  weUs 
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(Holmsen,  1959;  Rasmussen,  1974).  Mesozoic  volcanism  is  apparently  confined  to  S.  Scandia 
and  the  FFZ  (Bergstrpm,  1982;  Klingspor,  1976),  but  possibly  also  in  the  Egersund  Sub- 
Basin  (58N:6E)  (Fumes  et  al.,  1982).  Although  there  is  no  ultimate  evidence  of  Teniary 
volcanic  activity  in  Skagerrak,  there  are  indications.  For  example,  analysis  of  ash  layers  of 
Paleocene-Eocene  ages  in  Denmark,  imply  that  likely  source  areas  should  be  located  within 
100-200  km  bearing  N-NE  from  Limfjord  (LF  in  Fig.  8).  A  few  layers  also  indicate  a  NW- 
W  bearing  according  to  Nielsen  and  Heilmann-Clausen  (1988),  although  the  magmatic 
source(s)  have  not  been  found.  Other  investigations,  notably  using  geophysical  data  (Am, 
1973;  Hovland,  1987)  and  dredging  (Noe-Nygaard,  1967)  give  candidate  locations  as  marked 
on  Fig.  8  with  age  ranges  from  Permian  to  Tertiary. 

In  summary,  evidence  of  volcanic  activity  throughout  Carboniferous  to  present  have 
been  found  and  is  essentially  confined  to  areas  affected  by  rifting  and  subsidence,  and  hence 
to  areas  with  thin  crust  (thicknesses  less  than  about  30  km).  West  coast  of  Sweden  and  the 
northern  part  of  OG  (Fig.  8)  is  exeptional  here  with  crustal  thicknesses  of  approximately  36 
km. 


5.2  CRUSTAL  VELOCITY  -  DEPTH  DISTRIBUTIONS 

Crustal  velocity  plots  from  refraction  and  wide-angle  lines  (Fig.  9)  show  first  and/or 
second  order  discontinuities  in  central  parts  of  the  crystalline  crust.  On  the  Larvik  -  Lysekil 
profile  and  the  northern  part  of  E-4  there  are  two  discontinuities,  the  upper  one  at  approx. 
14  and  12km  respectively.  Even  though  the  P-velocity  depth  distribution  in  the  crust  is  not 
well  constrained  in  general  this  depth  range  coincide  with  the  onset  of  laminae  (Larsson  and 
Husebye,  this  volume),  tentatively  taken  to  mark  the  transition  from  the  upper  brittle  to  the 
lower  ductile  crust.  Similar  results  have  been  reported  by  Wever  (1989).  It  seems  then  that 
the  upper  and  lower  crust  is  roughly  equal  in  thickness  here.  Relativ  pronounced  velocity 
discontinuities  are  found  along  the  Fennolora  and  the  Fedje  -  Grimstad  profiles.  Areas  with 
intermediate  velocities  (7.3-7. 9  km/s),  that  are  higher  than  normal  crustal  velocities  and  less 
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than  typical  sub-Moho  velocities,  are  located  in  areas  with  truly  thick  crust  along  the 
Fennolora  profile.  Similar  features  have  been  observed  for  SW.  Norway  (S.  Deemer, 
pers.comm.). 

53  EARTHQUAKE  (EO)  LOCATIONS  -  CRUSTAL  THICKNESSES 

The  Nordic  EQ  database  (Seismological  Observatory,  Helsinki)  has  been  used.  We 
have  briefly  examined  the  epicenter  distribution  versus  our  estimated  crustal  thicknesses.  The 
time  period  covered  by  this  database  is  from  1375  to  present,  a  natural  subdivision  here  is: 
i)  Historical  data:  1375  -  1889.  ii)  Macro  seismic  data:  1890  -  1950.  iii)  Instrumental  data: 
1951  -  present.  For  example,  historical  data  are  generally  of  poor  quality  with  most  occuring 
EQ  left  unreported  and  at  least  larger  EQ  having  unreliable  epicenter  locations.  Macroseismic 
data  are  of  better  quality  reflecting  information  on  felt  EQ  through  circulation  of 
questionaries  (e.g.,  see  Husebye  et  al.,  1978).  Instrumental  data  are  relative  superior  to  the 
other  data  sets,  particularly  regarding  detectability  of  small  events,  and  also  those  located 
off-shore.  Both  historical  and  macroseismic  epicenter  locations  are  landward  biased.  Among 
the  four  hypocenter  parameters;  focal  depth  is  the  least  reliable,  and  often  a  default  value  of 
15  km  is  used  even  for  instrumental  data. 

This  being  said,  we  have  plotted  epicenter  locations  as  a  function  of  crustal  thickness 
both  on  a  map  (Fig.  10)  and  also  in  a  histogram  fashion  (Fig.  11).  Indeed,  there  is  a  striking 
correlation  between  EQ  locations  and  crustal  thicknesses  which  are  even  more  pronounced 
for  large  magnitude  events.  A  high  EQ  concentration  in  W.  Norway  and  off  the  west  coast 
is  obvious  from  Fig.  10.  There  are  also  some  events  in  south-central  Sweden  but  these  are 
mainly  of  small  magnitude  (Ml  <  3.0).  Most  events  are  located  in  areas  relatively  recently 
subjected  to  extension  where  the  crust  is  thin  and/or  the  crustal  thickness  gradient  is  high, 
this  is  even  more  striking  for  the  largest  earthquakes.  The  same  phenomenon  is  observed  for 
continental  crust  on  a  global  scale  as  demonstrated  by  Johnston  (1989). 
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5.4  CRUSTAL  THICKNESS  -  TECTONIC  AGE 

The  crustal  thickness  appears  to  be  increasing  as  a  function  of  tectonic  age.  The 
crustal  thickness  has  been  sampled  evenly  and  using  the  age  tied  to  the  latest  period  of 
major  tectonic  activity  for  a  given  area.  For  the  Archean  period  the  data  have  been  taken 
from  Luosto  and  Korhonen  (1986),  giving  rather  few  points  compared  to  the  other 
periodes/areas.  Fig.  13  and  Table  V  contain  the  results  for  the  area  investigated  and  reveal  a 
marked  thickening  of  the  older  parts  of  the  shield,  except  for  a  minor  thinning  in  the 
(oldest)  Archean  province.  This  has  been  observed  in  Australia  as  well  (Drummond  and 
Collins,  1986).  Meissner  (1986)  has  expressed  crustal  thickness  as  a  function  of  age  using 
the  relation:  2^251ogt-32,  where  Z=thickness  and  t=age  in  Ma,  using  global  data.  A  linear 
function  like  Z=8i+31  where  Z=thickness  and  t=age  in  Ga  gives  a  reasonable  fit  to  our  data, 
as  shown  in  Fig.  13  together  with  that  of  Meissner.  It  appears  that  after  about  1500  Ma  the 
continental  growth  terminated,  that  is  accretion  of  microcontinents  and/or  continental 
fragments  essentially  ceased  (e.g.,.  see  Berthelsen  and  Marker,  1986).  There  also  appears  to 
be  a  limit  to  crustal  thickness  growth  over  time  as  discussed  by  Warren  (1989).  The  essence 
of  current  hypothesis  is  that  during  crustal  cooling  in  shield  areas  an  instability  is  created  in 
the  lower  crust  in  terms  of  a  modest  density  increase  due  to  phase  transformation.  The  effect 
is  a  presumed  mass  exchange  between  the  lower  crust  and  the  upper  mantle,  thus  limiting 
crustal  growth  to  a  thickness  of  50  -  55  km. 

5.5  GRAVITY 

The  long  wavelength  Bouguer  gravity  map  of  the  area  is  displayed  in  Fig.  14 
(Balling,  1984).  It  exhibit  a  generally  smooth  picture  with  some  positive  feature  in  northern 
Skagerrak  and  one  of  the  Swedish  west  coast  with  a  NW  trend,  see  also  Wolf,  1984.  The 
dominating  feature  is  a  marked  negative  anomaly  (  -80  mgal)  over  western  Norway 
Caledonides,  and  thus  coincides  with  sub-Moho  Pn  and  Sn  low  velocity  anomalies  (Bannister 
et  aL,  this  volume).  This  and  the  fact  that  the  gravity/crustal  thickness  are  poorly  correlated 
in  this  area  inplies  a  sub-Moho  mass  deficit.  The  origin  of  this  mass  deficite  is  tentatively 
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tied  to  the  Caledonide  collision  per  se  and/or  the  opening  of  the  Norwegian  Sea  58  Ma  ago. 

6.  SUMMARY  AND  CONCLUDING  REMARKS 

The  oldest  part  of  the  Baltic  Shield,  east  of  the  PZ,  is  thick  and  relatively 
undisturbed.  The  Moho  depth  and  the  crustal  thickness  varies  from  40  to  50  km,  and  their 
gradients  are  rather  small,  with  the  exception  of  the  southernmost  part  of  Sweden. 

In  the  southern  part  of  Kattegat  the  FFZ/FBZ  plays  an  important  role,  and  these 
trends  are  dominating  both  the  Moho  depth  and  the  crustal  thickness  map.  The  Moho  depth 
varies  from  32  to  36km  across  the  zone  and  its  gradients  are  steep.  This  is  even  more 
pronounced  on  the  crustal  thickness  map.  The  FFZ/FBZ  looses  importance  westwards  as  the 
Moho  depth  and  crustal  thickness  variations  get  smaller,  and  the  FFZ/FBZ  covers  a  broader 
area. 

In  south  Skagerrak  the  NDB  seems  to  be  the  dominating  feature  with  a  generally 
elevated  Moho,  thin  crust,  small  gradients  and  no  dramatic  changes  in  the  westward  direction 
before  the  CG/VG.  However,  the  subsidence  history  north  and  south  of  the  FFZ  is  different 
as  demonstrated  by  Pedersen  et  al.  (this  volume). 

Further  north  in  Skagerrak  the  Moho  depth  and  crustal  thickness  are  dominated  by 
the  OR/OG.  The  Moho  elevation  is  about  5  km.  The  OR  appears  to  be  divided  in  two  parts, 
a  northern  dominated  by  volcanic  activity  and  large  amounts  of  volcanic  material  extruded 
during  a  rather  short  period  of  40-60  Ma  (Neumann  et  al.,1986).  The  large  volumes  of 
extrusives  and  the  elevated  Moho  would  indicate  a  mantle  source  as  the  prime  driving  force 
for  the  OR.  This  also  agrees  with  the  models  for  mantle  plumes  put  forward  by  Neugebauer 
(1983).  Also  Lynch  and  Morgan  (1987)  argues  that  the  prime  factor  controlling  rifting  of 
continental  crust  is  variation  in  the  geotherm. 

The  lateral  movements  along  the  FBZ/FFZ  is  estimated  by  several  authors  to  be  in 
the  order  of  10-20  km  (Liboriussen  et  al.,1987;  Anonymous,  1988).  These  movements  may 
have  controlled  the  southern  part  of  the  OR  where  large  scale  block  faulting  have  been 
observed  (Husebye  et  al.,  1988);  see  also  Fig.  4  and  5. 
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The  tounge  of  thicker  crust  at  the  west  coast  of  Norway  coincides  with  the 
Caledonides.  It  terminates  towards  the  CGA^G  in  the  North  Sea.  A  similar  feature  is  found 
on  the  British  side  of  the  CG/VG  as  mentioned  above. 

Volcanic  activity  throughout  Carboniferous  to  present  is  confined  to  areas  affected  by 
rifting  and  thin  crust. 

Velocity  discontinuities  coincides  with  the  onset  of  laminae,  taken  to  mark  the 
transition  between  the  upper  brittle  and  the  lower  ductile  crust. 

EQ  occurrences  show  that  most  events  are  located  where  the  crust  is  thin  and/or  the 
gradient  is  steep.  This  is  even  more  pronounced  for  the  high  magnitude  (ml>3.0)  events. 

Cmstal  thickness  increases  with  age  until  about  2000  Ma  and  then  level  off  or  is 
slightly  diminishing,  and  apparently  no  new  crust  has  been  accreted  after  about  1500  Ma. 
However,  crustal  thicknesses  seemingly  reflect  tectonic  age  as  the  Sveconorwegian  (1200  - 
850  Ma)  crust  is  relatively  thin  being  deformed  in  post  Cambrian  times. 

The  long  wavelength  Bouguer  gravity  have  a  marked  negative  anomaly,  this  coinsides 
with  the  Caledonides.  Together  with  Pn  and  Sn  low  velocity  anomalies  this  indicate  a  sub- 
Moho  mass  deficit. 

This  study  has  been  aimed  at  presenting  an  overview  over  currently  available 
geological  and  geophysical  (seismic)  information  bearing  on  the  crustal  structures  in  S. 
Scandinavia  and  adjacent  areas.  In  particular,  crustal  thickness  and  Moho  depths  over  this 
region  have  been  contoured  at  the  2  km  level.  There  is  an  obvious  correlation  between  the 
various  geological  and  geophysical  features  but  on  the  other  hand  such  relationships  are  most 
difficult  to  quantify.  We  have  not  attempted  to  do  this  for  the  simple  reason  that  the 
underlying  geodynamic  processes  deforming  the  crust/lithosphere  are  not  well  understood. 
These  are  indeed  intriguing  problems  where  progresses  have  been  slow  despite  much 
research  efforts  over  the  last  decade.  Hence  we  limited  the  scope  of  our  study  to  presenting 
present  crustal  information  which  we  consider  to  be  a  good  starting  platform  for  geodynamic 
modelling  exercise  now  being  contemplated. 
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TABLE  I 


TECTONIC  EVOLUTION  OF  THE  BALTIC  SHIELD 
-  MAJORE  DEFORMATION  PHASES 


OROGEN/AREA 

TIME  (MA) 

MAIN  FEATURES 

Saamian 

3100-2900 

Basement  of  late  archean  greenstone  belt. 

Lopian 

2900-2650 

Generated  large  volumes  of  new  crust. 

Svecofennian 

2000-1750 

Large  volumes  of  igneous  rocks  (1880  ±20  Ma). 

Trans  Scand. 
Igneous  Belt 

1785-1650 

Ensialic  spreading. 

Sveconorwegian 

Gothian 

1750-1500 

Hot  lithospheric  regime,  building 
the  SW.  Scandinavian  region. 

Hollandian 

1500-1400 

"Interlude"  with  mafic  intrusions. 

Sveconorwegian 

Greenvillian 

1250-  900 

Rotation  of  the  Shield,  complex  orogen 
with  large  scale  reworking  of  old  crust. 

Caledonian 

600-  400 

Collision,  trust-Nappe  tectonics. 

Oslo  Graben 

Oslo  Rift 

305-  245 

Rifting,  large  volumes  of  volcanic 
material,  large  normal  faults. 
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TABLE  11 


EVOLUTION  OF  DENMARK  AND  ADJACENT  SEAS 


PERIOD  MAIN  FEATURES 

Precambrian  Few  details  due  to  lack  of  sediments. 

Early  Paleozoic  Thick  sediment  strata  in  grabens  to  the  NE. 

Silurian  Basin  subsidence  to  the  SW  of  FBZ. 

Devonian  Most  of  Denmark  uplifted  -  erosion.  Strong  faulting  and 

-  Carboniferous  uplift  of  marginal  basins.  The  Caledonian  orogen  of 

-  Hercynian  moderate  influence  -  Lower  paleozoic  strata  show  no 
evidence  of  compressive  deformation. 

Late  Carboniferous  FBZ  activated,  dextral  movement,  wrench  tectonics; 

-  Early  Permian  complex  faulting  and  widespread  magmatic  activity. 

Rice  of  RFH  and  erosion.  Taphrogenesis  -  rifting 
initiated  in  Rpnne,  Oslo  and  Horn  Grabens. 

Early  Permian  Volcanic  activity  abated.  NDB  subsidence. 

RFH  a  positive  structure. 

Triassic  Repeated  reactivation  of  FBZ  -  transtensional  regime. 

Denmark  -  crustal  extension,  rapid  differential  subsidence, 
decreasing  in  late  Triassic.  Block  faulting  and  subsidence 
along  FBZ.  Mantle  derived  basaltic  plugs  onshore  Scandia. 

Cretaceous  FBZ  reactivated  several  times.  Volcanic  activity. 

Accumulation  of  thick  sediments  in  Alborg  and  R0nne  grabens. 

Late  Cretaceous  FBZ  inversion  with  a  shift  from  Early  transtensional 

-  Tertiary  to  transpressional  regime.  Partial  uplift  often 

related  to  the  Alpine  orogen. 
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TABLE  III 


VELOCITIES  USED  FOR  DEPTH  CONVERSION. 


Velocity  (km/s) 

Layer, T’eriod 

Lilhology/comment 

1.46 

Water 

Brackish,  Cold 

1.8 

Quaternary 

Clay,  Silt,  Sand,  Gravel 

3.8 

Upper  Cretaceous 

Chalk 

3.0 

L.Cretaceous-U.Jurassic 

Shale/Sand 

3.9 

L.Jurassic,  Triassic 

Shale/Sand 

5.3 

Upper  Permian 

Evaporites 

5.05 

Paleozoic 

Sediments 

6.1-7.0 

Precambrian 

Crystalline  crust 

8.1 

Sub  Moho 

The  velocities  are  based  on:  Wells:  Skagerrak-1  (57  50.6N  9  55E),  Skagerrak-l 
(57  26N  8  33. IE)  and  North  Sea-3  (57  0.8N  6  lO.lE);  refraction  experiment  from 
the  area  (Egilson  and  Husebye,  1989);  Faleide,  (1984)  and  processing  results  from 
the  M/V  Mobil  Search  data  (J.E.  Lie  pers.comm.). 
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TABLE  IV 


SEISMIC  PROFILING  -  CRUSTAL  MAPPING  STUDIES 


Map- ref 

Name/Area 

Refrence 

Type 

A 

Central  Graben 

Banon  and  Wood,  1984 

Refrac 

B 

Mobile  Search/ 

West  coast  of  Norway 

S.  Deemer  pers,  comm. 

W.Refl/ 

Refrac 

C 

M0re  margin 

Olafsson,  1988 

ESP 

D 

Cannobe/S.  Norway 

Cassell  et  al.,  1983 

Refrac 

E1-E5 

Eugeno-S/  Denmark, 
Kattegat,  SW. Sweden 

Anonymous,  1988/ 

Lund  et  al.,  1987 

Refrac 

F1-F3 

Fennolora  S.  pan 

Clowes  et  al.,  1987 

Refrac 

G 

Swedish  Lappland 

BSth,  1984 

Refrac 

H 

T  rondheim-S  unds  vail 

Vogel  and  Lund,  1971 

Refrac 

I 

Oslo  Trondheim 

Kanestr0m,  1971 

Refrac 

J 

Otta  -  Arsund 

Mykkelstveit,  1980 

Refrac 

K 

Flora  -  Asnes 

Sellevoll  and  Warrik,  1971 

Refrac 

L 

Fedje  -  Grimstad 

Sellevoll  and  Warrik,  1971 

Refrac 

M-N 

Oslo  Rift 

Tryti  and  Sellevoll,  1977 

Refrac 

O 

Larvik  -  Lysekil 

Egilson  and  Husebye,  1989 

Refrac 

P1-P4 

NSDP84-01:04 

North  Sea 

Fichler-Fettig  and 

Hospers,  1988 

Reflec 

NOA 

NORSAR,  Norway 

Berteussen,  1977 

Spec.r. 

COP 

HFS 

KON 

UME 

UPP 

Copenhagen,  Denmark 
Hagfors,  Sweden 
Kongsberg,  Norway 
Ume3,  Sweden 

Uppsala,  Sweden 

Bungum  et  al.,  1980 

Spec.r. 

Siljan,  Sweden 

Lund  et  al.,  1988 

Reflec 

S .  S  weden-Kattegat 

Komfajit  and  Larsson,  1987 

Reflec 

Viking  Graben 

Hospers  and  Ediriweera,1988  Reflec 

(Reftac  =  refraction  profiling;  W.Refl  =  wide  angle  reflection  profiling;  Spec.r. 
long  period  seismic  spectral  ratio  technique;  ESP  =  cxspanding  spread  profile.) 
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TABLE  V 


CRUSTAL  THICKNESSES  VS  AGE. 


AREA 

TIME  (MA) 

AVERAGE 
THICKNESS  (km) 

ARCHEAN 

2700-2500 

49.3 

SVECOKARELIAN 

2100-1900 

51.5 

SVECOFENNIAN 

2000-1800 

46.0 

TRANS  SCAN.  IGNEOUS  BELT 

1750-1550 

44.3 

SVECONORWEGIAN 

1200-  850 

37.6 

SN  EAST  of  OG 

39.7 

SN  WEST  of  OG 

34.5 

CALEDONIAN 

550-  400 

37.3 

CALEDONIAN 

38.6 

WESTERN  GNEISS  R. 

34.1 

OSLO  GRABEN 

305-  245 

31.5 
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Figure  Captions 


Fig.  1  Location  map  for  the  Skagerrak  profiles  shot  by  MA^  Mobile  Search  in  the  winter 
1987  (Husebye  et  al.,1988). 


Fig.  2  Structural  map  of  south  Scandinavia  showing  main  features.  CF:  Caledonian  Front; 
CG:  Central  Graben;  ESB:  Egersund  Sub-Basin;  FBZ:  Fennoscandian  Border  Zone;  FSE: 
Fasrpy-Shetland  Escarpment;  FFZ:  Fjerritslev  Fault  Zone;  HG:  Horn  Graben;  KBFC: 
Kristiansund-Bod0  Fault  Complex;  MB:  M0re  Basin;  MTFZ:  M0re-Tr0ndelag  Fault  Zone; 
MZ;  Mylonite  Tjont]  OG:  Oslo  Graben;  PZ:  Protogine  Zone;  RFH:  Ringk0bing-Fyn 
High;  SF:  Svecofennides;  SNF;  Sveconorwegian  Front;  TEF:  Trans  European  Fault;  VG: 
Viking  Graben;  WG:  Western  Gneiss  Region.  (Modified  after  EUGENO-S  Working  Group, 
1988;  Bukovics  and  Ziegler,  1985;  Gorbatschev,  1985;  Rawson  and  Riley,  1982).  (These 
abreviations  are  frequently  used  in  the  text). 


Fig.  3  Seismic  profiles  and  seismograph  stations  in  S.  Scandinavia  used  as  sources  for 
crustal  thickness  estimates  (profiling  details  and  references  in  Table  IV). 


Fig.  4  Linedrawings  of  OG-2  and  OG-12.  Depth  in  km.  Intersecting  lines  are  marked. 


Fig.  5  Linedrawings  of  OG-13,  OG-5  and  OG-7.  Depth  in  km.  Interseaing  lines  are 
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marked. 


Fig.  6  Moho  depth  below  mean  sea  level  -  2  km  contour  interval. 


Fig.  7  Thickness  of  crystalline  crust  in  km  -  2  km  contour  interval. 


Fig.  8  Volcanic  activity  and  thickness  of  crystalline  crust  in  km  -contour  interval  4  km. 
Arrows  indicate  directions  to  early  Tertiary  magmatic  sources  within  150  km  from  ash 
deposits  at  Limfjorden  (LF),  lines  are  Permian  dykes,  circles  are  Permian  volcanism, 
diamonds  are  Mesozoic  volcanism  and  stars  are  possible  Tertiary  volcanism. 


Fig.  9  Velocity  plots  a:  Larvik  -  Lysekil  and  northen  part  of  EUGENO-S-4:  b:  Fedje  - 
Grimstad;  and  c:  Fennolora,  southern  part. 


Fig.  10  Epicenters  in  S.  Scandinavia  1951  -  1988,  Open  squares:  Ml;2.4-3.5;  filled 
squares:  Ml:3. 6-5.5  and  crustal  thickness  in  km,  contour  interval  2  km. 

Fig,  11  EQ  frequencies  vs  crustal  thickness.  Open:  Ml:2.4-3.5  ;  Hatched:  Ml:3.6-5.5. 

Fig.  12  EQ  frequencies  vs  crustal  thickness,  adjusted  for  unit  area  of  crustal  thickness 

(approx.  25x25km^).  Open:  Ml:2.4-3.5  ;  Hatched:  Ml;3.6-5.5. 
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Fig.  13  Crustal  thickness  vs  tectonic  age.  Stars  =  average  for  period;  Triangles 
average  for  subgroup.  See  Table  V  for  details. 

Fig.  14  Bouguer  gravity  map,  contour  interval  20  mgal  (Balling,  1984). 
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Abstract 

Pn  and  Sn  velocity  estimates  are  derived  for  sub-Moho  structure  beneath  Fennoscandia 
using  a  tomographic  conjugate  gradient  scheme.  Observational  data  stem  from  local 
earthquake  recordings  of  crustal  phases  Pg,  Sg  (Lg)  and  subcrustal  phases  Pn  and  Sn  by 
the  Fennoscandian  seismograph  network.  Unknowns  arc  Pn  and  Sn  velocities  along  with 
event  hypocentral  parameters.  Crustal  thicknesses  arc  presumed  known  and  arc  not 
estimate  while  velocities  arc  allowed  to  vary.  The  most  prominent  P  and  S  velocity 
anomalies  arc  found  in  the  northwest  (in-land  Lofoten)  and  southwest  (M0rc)  of  Norway 
and,  for  shear  waves,  also  in  central  Hnland.  There  is  a  good  correlation  between  the 
velocity  anomalies  in  Norway  and  observed  negative  Bouguer  anomalies,  indicating  the 
presence  of  low-density  sub-Moho  material  in  these  areas.  Geologically  the  Pn  and  Sn 
velocity  anomalies  appear  to  be  associated  with  geodynamic  processes  tied  to  the  opening 
of  the  Norwegian-Greenland  Sea. 

1.  introduction. 

Seismology  currently  offers  the  most  powerful  means  for  mapping  the  physical 
properties  of  the  earth,  leading  towards  the  understanding  of  present  and  past 
geodynamical  processes.  Crustal  velocity  structures  have  traditionally  been  mapped  only 
along  selected  profile  lines  using  refraction  and  wide-angle  reflection  data  (e.g.,  in 
projects  such  as  FENNOLORA,  Galson  &  Mueller,  1986).  Although  such  experiments 
give  detailed  information  on  vertical  velocity  variations  the  resolution  of  lateral 
heterogeneities  is  restricted  to  the  profile  line.  Present  knowledge  of  3-D  velocity 
structures  therefore  comes  from  seismological  research  where  earthquake  signals  arc 
continuously  recorded  by  seismic  stations.  The  first  effort  in  this  field  was  that  of  Aid  et 
al.  (1977)  using  arrival  times  of  teleseismic  P-waves  (distances  >  3000  km)  at  NORSAR, 
Norway.  The  method  has  since  been  used  for  mapping  larger  scale  velocity  anomalies  in 
Fennoscandia  (Husebye  and  Hovland,  1982)  and  in  southern  Scandinavia  (Husebye  et  al., 
1986).  Dispersion  analysis  of  Rayleigh  waves  (e.g.  Knopoff,1972  and  Der  and  Landisman, 
1972)  also  provides  information  on  lateral  velocity  variation,  resolving  the  variation  of 
shear  wave  velocity  with  depth  and  hence  the  lithosphere  thickness.  Calcagnile  (1982) 
applied  this  method  to  the  Fennoscandian  region  and  found  large  thicknesses  (>170  Ion) 
in  Finland  and  northern  Sweden  while  normal  thicknesses  (ca.llO  km)  were  found  for 
Denmark  and  south-western  Norway  (Dost,  1986). 

In  this  paper  we  use  arrival  times  of  Pn  and  Sn  crustal  phases,  as  reported  by 
Fennoscan^an  stations,  to  obtain  tomographic  estimates  of  the  sub-Moho  P  and  S 
velocities  in  the  region.  Similar  studies  have  been  made  by  Aki  and  Lee  (1976)  and 
Hearn  and  Clayton  (1986a,b)  using  data  from  the  dense  Southern  Californian  seismic 
network.  Our  method  differs  from  that  used  by  Hearn  and  Clayton,  which  is  a  delay  time 
method,  in  that  we  allow  for  variations  in  crustal  thickness  and  include  the  earthqu^e 
source  parameters  as  unknowns  during  inversion,  thus  effectively  relocating  the  events. 
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We  also  use  a  different  numerical  algorithm,  results  are  presented  here  for  a  least  square 
conjugate  gradient  (CGLS)  algorithm  (e.g.  see  Spakman  and  Nolet,  1988). 

2.  Model  description  and  inversion  technique 

ITie  model  used  in  the  study  consists  of  a  two  layered  lithosphere  which  is  divided  by 
the  Moho  discontinuity,  the  upper  layer  containing  the  whole  crystalline  crust,  the 
critically  refracted  Pn  and  Sn  rays  travelling  through  the  lower  layer  representing  the  sub- 
Moho  region.  A  similar  two  layer  model  has  earlier  proved  sufficient  for  locating  local 
events  in  Fennoscandia  (Ruud  et  al.,  1988).  A  zero  velocity  gradient  was  assumed  in  the 
sub-Moho  region;  in  most  earth  models  very  low  or  even  negative  gradients  are  given  for 
the  earths  upper  100  km.  In  the  model  of  Herrin  (1968)  for  example,  which  is  a 
smoothed  model  with  no  discontinuities  or  low  velfv'ity  zones,  the  P  ray-path  lies  within 
the  upper  100  km  for  distances  of  up  to  ca  1500  km. 

Both  layers  were  divided  into  36  x  42  rectangular  blocks  of  constant  velocity,  each 
block  l°(long)  X  0.5°(Iat).  The  large  lateral  extent  of  the  model  requires  that  the  earth 
sphericity  is  taken  into  account.  The  block  subdivision  of  our  model  was  along  lines  of 
constant  latitude  and  longitude;  rays  which  are  assumed  to  follow  great  circles  will 
appear  slightly  curved  relative  to  the  model  grid  system.  All  lateral  distance  calculations 
were  performed  in  spherical  coordinates. 

Variations  in  crustal  thickness  were  allowed  for  by  varying  the  thickness  of  the  upper 
layer  blocks  (but  everywhere  with  horizontal  interface  to  the  block  below  in  the  lower 
layer).  Detailed  knowledge  of  the  variation  of  the  Moho  depth  beneath  the  region  has 
been -obtained  from  deep  reflection  seismic  work;  it  is  known  that  the  Moho  depth  varies 
by  more  than  30  km,  from  ca.  20  km  off  the  west  coast  of  Norway  to  some  50  km  in 
central  Finland.  The  a-priori  model  of  crustal  thickness  used  in  this  study  (Fig  1)  was 
based  on  the  works  of  Kinck  et  al.(1989,  this  volume),  Luosto  (1989,  this  volume), 

Planke  et  al.(1989,  this  volume),  Faleide  et  a].(1989)  and  Prodehl  and  Kaminski  (1984). 

Ray-tracing  was  carried  out  initially  using  a  model  with  no  lateral  velocity  variation. 

The  rays  were  traced  as  straight  lines  between  the  hypocenter  and  station  down  to  the 
Moho  interface,  using  the  angle  of  incidence  calculate  from  the  reference  layer  velocities 
(while  ignoring  the  effect  of  potential  slanting  Moho).  Between  the  two  interface 
intersection  points  ray-paths  were  assumed  to  be  horizontal  and  completely  within  the 
lower  layer.  Differences  in  the  level  of  the  two  intersection  points  were  accounted  for  by 
applying  a  length  scaling  factor  to  the  segment  of  the  ray  in  the  lower  layer.  P-wave 
velocities  of  6.32  and  8.18  km/s  were  used  as  initial  reference  velocities  for  the  crust  and 
sub-Moho  respectively;  initial  S  velocities  were  3.65  and  4.65  km/s.  These  reference 
velocities  were  deriv^  from  a  least  squares  regression  analysis  of  our  data  set  described 
below. 

2.1  Inversion 

Velocity  anomalies  and  source  parameters  are  highly  correlated  in  travel  time 
tomography  involving  sources  which  are  not  precisely  located.  Both  types  of  variable 
should  be  treated  as  unknowns  during  the  estimation  process  (Aki  and  Lee,  1976; 

Spencer  and  Gubbins,  1980).  We  thus  considered  the  travel  time  residuals  to  consist  of 
three  parts;  (1)  that  caused  by  genuine  velocity  anomalies,  (2)  mislocation  of  the  events, 
and  (3)  errors.  The  latter  group  consists  of  not  only  reading  and  timing  errors  but  also 
errors  resulting  from  deficiencies  of  the  simplified  model  such  as  the  artificial  block 
structure  and  the  lack  of  vertical  velocity  gradient. 


The  system  of  equations  can  then  be  written,  following  Spakman  and  Nolet  (1988), 


Ax  =  [L,  G]x  =  d  +  e, 


(1) 


where  s  is  the  vector  of  slowness  anomalies,  g  the  vector  of  hypocentre  parameter 
corrections  (origin  time,  latitude,  longitude  and  depth),  d  the  vector  of  time  residuals,  e 
the  error  vector,  G  the  matrix  of  partial  derivatives  with  respect  to  hypoccntral 
parameters,  and  L  the  matrix  of  ray  segment  lengths.  Station  corrections  are  accounted 
for  within  the  upper  layer  block  on  which  each  station  is  situated.  Assuming  that  the 
errors  are  random  and  approximately  follow  a  gaussian  distribution  with  zero  mean  we 
want  to  obtain  a  least  square  estimate  of  the  x  vector.  Appropriate  damping  and 
weighting  (according  to  a-priori  information)  was  applied  to  the  system  (1)  as  detailed  by 
Spakman  and  Nolet  (1988).  Such  weighting  is  essential  for  a  meaningful  result  as  the 
vector  X  of  unknowns  in  (1)  consists  of  different  physical  parameters. 

In  our  study  the  system  (1)  consists  of  thousands  of  unknowns  and  equations  so  that 
solution  by  generalized  matrix  inversion  is  prohibited  by  computer  time  and  memory 
considerations.  We  therefore  seek  a  solution  using  row-action  or  projection  techniques. 
These  are  iterative  algorithms  with  modest  memory  requirement  resulting  in  stable 
solutions  after  relatively  few  iterations.  One  such  ^gorithm  we  have  applied  is  the  least 
square  conjugate  gradient  (CGLS)  (Hestenes  and  Stiefel,  1952).  This  algorithm  produces 
results  very  similar  to  the  popular  LSQR  method  (Paige  and  Saunders,  1982).  We  also 
carried  out  experiments  using  a  SIRT  algorithm  (e.g.  Nolet,  198B).  This  later  algorithm 
tends  to  have  errors  in  poorly  sampled  blocks  and  was  not  used  for  the  results  presented 
below. 

As  tomographic  systems  arc  usually  under-determined  (at  least  in  some  parameters)  it  is 
critical  that  the  solution  is  stabilised.  Stabilisation  in  CGLS  results  from  the  use  of 
explicit  damping  together  with  implicit  damping  within  the  actual  algorithms,  the  smaller 
eigenvalues  gradually  entering  into  the  solution  with  increased  iterations  (Spakman  and 
Nolet,  1988).  A  mild  smoothness  constraint  was  also  imposed  during  the  reconstruction 
process  (e.g.  Herman  ct  al.,  1975)  using  a  nearest  neighbour  moving  average  filter  every 
fifth  iteration.  Such  smoothing  assists  convergence  towards  a  solution  involving  smoothly 
varying  anomalies. 

3.  Data  and  their  inversion 

Arrival  times  as  reported  from  seismological  stations  in  Fennoscandia  (Fig  2)  arc 
routinely  collected  by  the  Seismological  Institute  at  the  University  of  Helsinki,  Finland, 
where  event  association  and  epicenter  location  arc  also  carried  out.  Data  from  1981  to 
the  present  is  available  in  a  computer  readable  form  although  the  reporting  frequency  and 
station  coverage  were  somewhat  incomplete  prior  to  1984.  This  data  base  is  small 
compared  to  the  one  ttsed  by  Hearn  and  Claytuu  (i986a,b)  but  resolution  is  fairly  good 
in  the  central  part  of  the  region  where  cross-sampling  is  best,  as  demonstrated  below; 
there  is  a  reasonable  distribution  of  events  across  the  region,  with  a  high  density  in  the 
west  Norway  and  Finland  areas  as  well  as  in  the  Norwegian  Sea. 

We  used  Pn  and  Sn  arrival  times  reported  from  seismic  stations  (Fig  2)  in  the 
Fennoscandia  region  between  January  1984  and  May  1988.  The  network  at  hand  is  too 
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coarse  for  incorporatirg  Pg  and  Sg(Lg)  phase  observations  as  these  phases  arc  seldom 
reponed  reliably  beyond  epicenter  distances  of  2-300  km  due  to  Pn  and  Sn  interference. 
The  arrivals  were  extracted  by  subjectively  ’windowing’  on  a  travel-time  distance  plot 
(e.g.  Fig  3)  using  an  absolute  residual  time  threshold  of  3s;  the  rejected  arrivals  scattered 
outside  of  these  windows  are  likely  to  be  associated  with  events  that  arc  poorly  located 
or  involve  data  input  errors.  More  than  1800  Pn  and  1500  Sn  arrivals  were  extracted 
between  distances  from  350  to  1500  kms.  Fig  4  shows  the  spatial  coverage  of  the  Pn  and 
Sn  raypaths  respectively.  It  is  seen  that  there  is  a  high  density  of  raypaths  in  the  center 
of  Fennoscandia,  along  the  west  coast  of  Norway  and  between  north-west  Norway 
(offshore)  and  the  Finland  seismograph  network. 

All  events  reported  for  the  time  period  were  relocated  by  us  both  before  the  inversion, 
using  the  updated  crustal  model  which  incorporates  the  variation  in  depth  of  the  Moho 
across  the  region  (Fig  1),  and  after  the  inversion,  incorporating  the  lateral  variations  of 
velocity  detailed  below.  TTie  epicenters  of  the  events  before  the  first  relocation  (i.e.  as 
reported  in  the  Finland  data  base)  are  shown  in  Fig  5.  together  with  the  final  locations 
following  the  completion  of  the  inversion  procedure.  Maximum  epicenter  shifts  from  the 
initial  ’Finland’  event  locations  to  the  final  inversion  locations  were  of  the  order  of  1® 
for  events  at  the  pheriphery  of  the  model,  but  of  the  order  of  10  km  for  events  in  its 
central  parts.  Most  of  this  epicentral  shift  occurred  in  the  first  pre-inversion  relocation, 
the  shift  in  the  epicenters  following  this  first  relocation  (i.e.  to  the  final  epicenters 
determined  by  the  end  of  the  inversion  procedure)  was  much  less  -  of  the  order  of  5  km 
in  the  central  areas,  approximately  half  of  the  above.  Most  of  the  epicentral  shift  can  be 
attributed  to  our  improvement  of  the  data  set  together  with  a  small  difference  in  the 
reference  Moho  velocity.  Event  depths  generally  shallowed  slightly,  remaining  within  the 
upper  crust  with  few  exceptions,  all  events  above  the  Moho.  Only  arrivals  from  events 
within  the  region  were  included  in  the  inversion  and  events  with  less  than  5  phases  were 
disregarded.  Station  elevation  corrections  were  applied  to  both  Pn  and  Sn  arrivals. 

Inversion  of  the  data  was  accomplished  through  application  of  the  COLS  algorithm, 
estimating  ca,1600  hypocenter  corrections  and  ca.l3()0  velocity  parameters.  The  RMS  of 
the  travel-time  residues  was  reduced  from  1.5  seconds  to  0.6s  (P-velocities)  and  1.5s  to 
0.74s  (S-velocities)  after  39  iterations.  A  mild  nearest-  neighbour  moving  average  was 
applied  to  the  resulting  velocity  images  to  emphasize  the  longer  wavelength  features 
(Herman  et  al.,  1975).  This  filter  had  a  small,  if  negligible,  effect  on  the  convergence.  A 
more  important  parameter  was  the  relative  scaling  between  hypocentral  and  slowness 
corrections,  discussed  further  below. 

3.1  Resolution  and  error  analysis 

In  the  following  section  the  major  anomalies  found  in  the  tomographic  inversion  are 
presented  and  discussed.  First,  however,  we  must  exanune  the  resolution  associated  with 
the  data  set  Use  of  the  row-action  techniques  prohibits  us  from  calculating  formal 
resolution  and  standard  errors  for  the  estimated  parameters.  As  a  substitute  we  inverted 
several  synthetic  data  sets  constructed  using  the  rays  from  the  actual  data  set  and  the 
travel  times  calculated  from  artificial  slowness  anomalies. 

The  synthetic  models  are  shown  in  Fig  6,  along  with  the  P-velocity  images  derived 
after  application  of  the  CX3LS  algorithm.  Tlie  first  synthetic  model  (Fig  6a)  involved 
patterns  of  low  (7.9  km/s)  and  high  (8.5  km/s)  bodies  in  a  sub-Moho  with  a  background 
velocity  of  8.2  km/s.  No  anomalies  were  initially  placed  within  the  modelled  crustal 
layer,  which  was  assigned  an  initial  velocity  of  6.32  km/s.  Velocities  of  all  blocks  in 
both  the  sub-Moho  and  crustal  layers  were  included  as  unknowns  in  the  tomographic 
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reconstruction  (Fig  6b),  together  with  hypocentral  correction  parameters,  which  were 
initially  set  to  zero.  The  second  model  (Fig  6c  and  6d)  involved  east-west  oriented 
velocity  anomalies  (7.9  krn/s  and  8.5  km/s)  within  the  sub-Moho. 

Inversion  of  the  data  was  accomplished  through  application  of  the  CXjLS  algorithm. 
The  RMS  of  the  travel-time  residuals  for  the  first  data  set  was  reduced  from  0.63  s.  to 
0.12  s.  after  39  iterations.  The  artificial  anomalies  in  the  center  of  the  model  (Fig  6a)  are 
well  resolved,  both  spatially  and  in  amplitude  (Fig  6b).  These  regions  correspond  to 
central  Fennoscandia,  that  is  Bothnian  Bay,  southwest  Finland,  the  east  and  center  of 
Sweden,  and  central  Norway.  We  can  see  some  smearing  of  the  anomalies  in  the  NW-SE 
direction  in  the  northern  part  of  the  model,  an  effect  that  is  best  seen  in  the  first 
synthetic  model  (Fig  6b)  but  that  is  also  reflected  in  Fig  6d  where  broadening  of  the 
northern  row  anomalies  can  be  seen.  This  effect  likely  results  from  the  biased  (NW-SE) 
direction  of  the  raypaths  in  this  part  of  Fennoscandia,  involving  waves  from  events  off- 
coast  NW  Norway  recorded  at  seismograph  stations  in  central  Finland.  In  addition  to  this 
effect  we  can  see  broadening  of  the  row  anomalies  to  the  west  of  Norway  (Fig  6d);  the 
resolution  of  the  images  are  highly  dependent  on  the  cross-sampling  of  the  ray-paths  (Fig 
4).  The  amplitudes  of  the  reconstructed  anomalies  are,  on  average,  ca  60  %  of  the 
synthetic  anomalies.  In  this  respect  we  note  the  dependence  on  the  relative  scaling  of  the 
hypocenter  parameters  with  respect  to  the  velocity  perturbations.  The  unknown  event 
corrections  are  usually  not  well  constrained,  and  can  easily  absorb  part  of  the  delay  time 
which  may  otherwise  be  explained  by  velocity  heterogeneities.  In  a  joint  inversion  as 
carried  out  here  this  has  a  imping  effect  on  the  estimation  of  the  sub-Moho  velocity 
heterogeneities  -  the  anomalies  tend  to  be  reconstructed  with  less  amplitude,  and  thus 
represent  minimum  estimates.  This  effect  is  most  dominate  in  regions  which  have  poor 
ray  coverage,  biasing  the  event  parameters  but  down-weighting  the  velocity  perturbation. 
To  a  lesser  extent  this  same  trade-off  also  occurs  between  the  estimated  velocities  of 
blocks  within  the  crustal  layer  and  the  sub-Moho,  again  resulting  in  down-weighted 
estimates  for  the  sub-Moho  velocity  perturbations.  The  reconstruction  would  act  similarly 
with  real  data  involving  data  errors;  in  this  regard  there  is  a  subtle  balance  between 
reconstructing  the  full  magnitude  of  the  cell  slowness  anomalies  and  the  degree  of 
influence  of  the  data  errors  on  the  velocity  estimates.  This  subjective  scaling  is 
determined  on  the  basis  of  the  amplitude  behaviour  of  the  components  of  the  solution  - 
the  hypocentral  corrections  and  velocity  estimates. 

We  note  that  this  synthetic  modelling  will  not  separate  the  artifacts  in  the  solution  that 
may  result  from  mis-representation  of  the  earth  structure  or  our  parameterisation.  The 
possible  effect  of  such  errors  on  the  velocity  estimates  arc  discussed  further  below. 

4.  Results 

Results  of  the  tomographic  inversion  of  the  data  set  are  shown  in  figures  7  and  8  for 
the  P  and  S-arrival  data  sets  respectively.  The  general  trend  seen  in  the  P  velocity 
reconstruction  is  towards  high  velocities  in  the  oldest  and  thickest  part  of  the  Baltic 
shield  and  low  velocities  in  Norway  and  the  offshore  areas.  It  should  be  noted  that  areas 
which  are  not  sampled  (e.g.  Barents  Sea)  will  retain  the  reference  values  for  the  P  and  S 
velocities.  The  most  prominent  features  in  the  P-velocity  image  are  the  low  velocity 
anomalies  beneath  western  (60-62°)  and  central  (62-63°)  Norway,  the  Lxifoten  area  in 
the  NW  (67-69°N)  of  Norway  and  the  weaker  anomaly  beneath  south-central  Finland  (see 
Fig  7  for  details).  Offshore,  the  negative  anomaly  trending  WNW  from  S.Sweden,  across 
Jutland  and  towards  SW  Norway  coincides  with  the  Fennoscandian  border  zone.  Weak 
high  velocity  anomalies  are  found  beneath  eastern  Sweden,  Bothnian  Bay  and  central  NE 
Finland.  From  the  synthetics  in  Fig. 6  the  above  anomalous  areas  are  well  sampled,  that  is 
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the  ray  path  coverage  is  good  with  respect  to  path  orientation  and  density  and  thus  the 
resolution  should  be  good.  However  the  results  from  the  synthetic  modelling  indicate  that 
the  anomalies  in  the  Lofoten  area  are  likely  to  be  smeared  slightly  in  the  NW-SE 
direction  while  those  anomalies  far  off-shore  western  and  southern  Norway  are  of 
marginal  significance. 

Similar  features  are  also  seen  in  the  S-velocity  image,  together  with  an  additional 
anomaly  at  65°N  (SE  Sweden).  Both  the  W. Norway  anomalies  are  quite  strong  and  again 
are  spatially  well  resolved.  The  S  velocity  anomaly  beneath  south-west  Finland  is  strong, 
involving  a  2-3%  contrast  with  the  surrounding  lithosphere,  slightly  greater  than  the  P- 
wave  anomaly  observed  in  that  area. 

Crustal  P  and  S  velocities  estimated  in  the  inversion  deviated  from  the  background 
reference  velocity  only  in  the  blocks  directly  beneath  seismograph  stations  and  events. 
Estimated  velocities  varied  from  6.17  km/s  to  6.45  km/s  across  the  region,  quite  a 
respectable  deviation  from  the  initial  6.32  km/s  but  well  in  line  with  the  range  of 
(average)  crustal  velocities  that  might  be  expected.  These  estimates  however  are  not 
further  interpreted  here  as  the  estimates  act  to  some  degree  as  buffers  for  data  errors,  as 
noted  above  in  the  synthetic  modelling. 

There  may  well  be  artifacts  resulting  from  mis-representation  of  the  earth  structure  or 
our  parameterisation  that  were  not  illuminated  by  the  synthetic  modelling.  For  example  if 
the  crustal  model  used  in  this  study  is  grossly  inaccurate  then  event  locations  and  the 
calculation  of  the  reference  ray  paths  will  both  be  affected  to  some  extent  although  the 
resulting  bias  in  the  images  is  difficult  to  detect  and  separate.  Such  errors  in  the  crustal 
model  are  not  considered  likely  however  as  the  crustal  model  in  Fig  1  is  quite  heavily 
sampled  apart  from  the  periphery  of  the  area,  and  is  based  on  many  intersecting  profile 
lines.  It  is  recognised  however  that  there  is  some  varying  interpretations  for  the  crustal 
model  in  south-west  Finland  and  for  the  crustal  transition  along  the  western  coastline  of 
Norway.  To  test  the  effect  of  possible  errors  we  carried  out  several  inversions  involving 
crustal  models  of  varying  thiclmess  (+/-  2  to  5  km),  and  found  only  slight  changes  in  the 
main  features  of  the  images  discussed  above.  For  example  when  the  model  crustal 
thickness  in  the  area  bounded  by  longitudes  10  and  13  and  latitude  63  and  66  is 
increased  by  5  km,  which  is  an  extreme  change  in  the  defined  nxxiel,  the  change  in  the 
P  velocity  estimate  for  the  sub-Moho  in  that  region  is  only  ca.0.04  km/s.  The  possible 
variation  in  the  sub-Moho  velocity  depends  on  the  sampling  of  the  particular  crustal 
block  in  question,  for  example  crustal  blocks  directly  above  the  western  coastline  of 
Norway  (the  region  delimited  in  the  above  example)  are  only  very  lightly  sampled.  The 
method  is  expected  to  be  most  sensitive  to  errors  in  the  crustal  thickness  beneath  stations 
although  these  errors  may  also  be  absorbed,  to  a  degree,  by  the  crustal  velocity  estimate 
for  the  particular  block.  Errors  in  the  thickness  beneath  epicenters  on  the  other  hand  will 
mainly  affect  the  estimated  hypocentral  depth. 

Other  errors  may  result  from,  for  example,  the  parameterisation  of  our  model  involving 
a  zero  gradient  within  the  upper  mantle,  and  the  subsequent  modelling  of  the  data  as 
seismic  energy  propagating  as  head  waves,  rather  than  diving  waves  within  the  upper 
mantle.  The  evidence  for  any  positive  velocity  gradient  in  the  mantle  beneath 
Fennoscandia  is  limited  and  there  is  generally  a  lack  of  accurate  information  on  the 
variation  in  P  and  S  amplitude  vs  distance.  Moreover  the  P  and  S  arrival  time  plots  (Fig 
3)  are  quite  linear  for  distances  up  to  around  1200  km,  suggesting  no  large  gradient  in 
the  mantle  lid.  If,  however,  a  large  positive  velocity  gradient  did  exist  beneath  the  region 
then  the  errors  induced  in  our  results  will  vary  depending  on  the  distance  of  propagation 
under  consideration,  being  higher  for  larger  distances,  and  presumably  leading  to  high 
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estimates  for  the  velocities.  It  is  difficult  to  quantify  the  overall  effect  on  the  velocity 
images  from  such  an  error  although  we  do  not  believe  that  it  would  be  great  -  75%  of 
the  arrivals  involve  epicentral  distances  of  less  than  1000  km.  Longer  epicentral  distances 
are  usually  associated  with  the  events  offshore  NW  Norway  and  corresponding  raypaths 
show  no  obvious  correlation  with  high  velocity  anomalies. 

5.  Geological  and  Geophysical  Framework  of  Fennoscandia 

The  anomalous  Pn  and  Sn  velocity  distributions  appear,  as  noted  above,  to  be 
correlated  with  dominant  geophysic^  features  in  the  region.  Before  discussing  this  in 
more  detail  we  first  outline  the  geological  and  geophysical  setting  of  Fennoscandia. 

The  Baltic  Shield,  comprising  Norway,  Sweden,  Finland,  the  Kola  peninsula  and  USSR 
Karelia,  was  mainly  formed  or  accreted  during  four  orogenies  between  3100  to  1500  Ma. 
The  shield  is  roughly  zoned,  in  the  sense  that  the  oldest  crust  is  found  in  NE  Finland 
and  the  Kola  peninsula  and  the  youngest  in  south  and  eastern  Norway  (Fig  9).  Most  of 
the  Baltic  shield  was  formed  during  the  Lopian  (2900-  2600  Ma),  the  Svecofennian 
(2000-1750  Ma)  and  the  Gothian  (1750-  1500  Ma)  orogenies  although  the  oldest  crust 
stems  from  the  Saamian  orogeny  (3100-  2900  Ma).  The  Protogine  zone  (PZ  in  Fig  9)  in 
south-central  Sweden  was  formed  by  episodic  faulting  accompanied  by  magmatic 
intrusions.  Today  it  constitutes  a  geologically  shaip  boundary  between  highly 
metamorphic  SW  Scandinavia  and  the  well  preserved  Svecofennian  terraiiis  to  liie  east. 
The  later  Sveconorwegian  orogeny  resulted  in  thorough  reworking  of  the  crust  as  far  east 
as  the  Protogine  zone  (Gaal  and  Gorbatschev,  1987).  The  Caledonian  orogeny  (600-400 
Ma)  involved  closing  of  the  lapetus  ocean  (Late  Ordovician)  and  subsequent  Himalaya- 
type  collision  between  Laurentia-Greenland  and  Fennosarmatia-Baltica.  Large  areas  of  the 
precambrian  crust  of  W. Scandinavia  become  concealed  beneath  basement-involving  nappes 
(e.g.  Gee  and  Sturt,  1985  and  Hossack  and  Cooper,  1986),  with  progressive  extensive 
metamorphism  in  a  NW  direction.  The  extent  of  crustal  deformation  is  not  known  in 
detail. 

More  recent  deformational  processes  have  involved  extensional  deformation  of  the 
continental  margin  offshore  Norway  and  substantial  thinning  of  the  crust  and  lithosphere. 
The  most  recent  deformation  process  of  interest  to  our  study  is  tl.e  opening  of  the 
Norwegian-Greenland  Sea  which  commenced  ca.  57.5  Ma  ago  (Talwani  and  Eldholm, 
1977),  and  the  transition  from  continental  to  oceanic  crust  tied  to  the  Faeroe-Shetland  and 
Vpring  Plateau  escarpments  offshore  (Fig.9).  As  the  extension  of  continental  lithosphere 
is  accompanied  by  influx  of  hot  magmatic  material,  anomalous  heat  flow  is  intuitively 
expected  in  certain  parts  off-shore  W.Norway  and  possibly  other  parts  of  the  continental 
margin.  This  in  turn  would  correspondingly  affect  the  Pn  and,  in  particular,  Sn  velocity 
distributions.  Key  references  on  the  geological  evolution  and  structural  mapping  of  these 
areas  are  Ziegler  (1988),  Planke  ct  al.(1989)  and  Faleide  et  al.(1989);  see  also  Fig  2  and 
Fig  9. 

5.1  Geophysical  investigations  -  crust/lithosphere  structural  mapping 

The  principal  geophysical  means  for  crusi/lithosphere  mapping  are  seismic  refraction, 
wide-angle  and  reflection  profiling.  The  main  results  within  the  region  are  incorporated  in 
the  crustal  thickness  model  for  Fennoscandia  (Fig  1)  and  have,  amongst  others,  been 
summarized  by  Kinck  et  al.(1989)  and  Luosto  (1989).  Tomographic  mapping  of  the  upper 
mantle  has  been  performed  by  Husebye  et  al.(1986)  using  teleseismic  P-wave  residuals, 
although  only  long  wave  length  structural  anomalies  could  be  mapped  using  this  data.  For 
an  overview  of  other  geophysical  investigations  see  Husebye  (’982)  and  others  in  this 


75 


volume.  'I'lie  gravity  surveys  are  of  most  interest  as  the  long  wavelength  anomalies  may 
reflect  sub-Moho  heterogeneities.  Tlie  Bouguer  gravity  map  for  Fennoscandia  is  displayed 
in  Fig  10.  One  outstanding  feature  is  the  pronounced  gravity  low  along  the  mountain 
range  of  Norway  which  broadens  over  the  south  central  and  Lofoten  areas.  Theilen  and 
Meissner  (1979)  and  Meissner  (1979)  noted  the  lack  of  mountain  roots  beneath  the 
mountain  range.  Tliey  suggested  that  at  least  part  of  the  implied  mass  deficit  must  be 
attributed  to  mantle  sources  while  smaller  wavelength  features  in  the  gravity  field  reflect 
the  more  superficial  Caledonian  nappe  lithology  (Wolff,  1983  and  Haworth  et  al.,  1988). 
Finally,  areas  exhibiting  high  heat  flow  are  often  characterized  by  reladvely  low  seismic 
velocities,  as  demonstrated  by  Hovland  et  al  (1981)  for  the  Central  Europe  lithosphere. 
Unfortunately  few  precise  heat  flow  measurements  are  available  for  Fennoscandia,  and 
they  are  particularly  lacking  in  the  younger  and,  for  this  study,  critical  ’coastal’  areas  of 
Norway  where  heat  imprints  from  the  opening  of  the  Norwegian-  Greenland  Sea  might 
exist. 

6.  Discus.sion 

TTie  dominant  shield  structures  in  central  parts  of  Fennoscandia  exhibit  only  modest  P 
and  shear  velocity  anomalies  (Fig  7  and  8).  This  is  not  unexpected  as  shield  areas  in 
general  are  seismologically  homogeneous.  The  notable  exceptions,  as  noted  above,  are  the 
S  and  P  velocity  anomalies  beneath  south  central  Finland  and  the  NE  coast  of  Sweden 
which  coincide  neatly  with  the  Svecofennides  (Fig  9).  Ihe  shear  velocity  anomaly  is 
dominant,  possibly  indicating  lateral  variarions  in  Poissons  ratio. 

The  other  dominant  velocity  anomalies  arc  those  found  in  Norway  which  coincide  with 
pronounced  negative  gravity  anomalies  of  long  wavelength  (Fig  9).  As  noted  above  the 
implied  mass  deficit  cannot  be  fully  explained  by  crustal  structure,  the  presence  of  low- 
density  upper  mantle  is  required  to  explain  the  deficit.  The  velocity  anomaly  determined 
for  the  Lofoten  area,  ca.  3%,  corresponds  approximately  to  a  density  anomaly  of  0.1-  0.2 
Mg/m^  (using  the  reladon  of  Birch,  1961).  A  body  10  x  200  x  200  km^  and  with  a 
relauve  density  contrast  of  O.lMg/m^  would  give  rise  to  a  Bouguer  anomaly  of  30-40 
mgal  or  roughly  half  of  the  observed  anomaly.  It  is  puzzling  however  that  critical 
sections  of  the  prolonged  negative  Bouguer  anomaly  beneath  the  Caledonides  do  not  have 
an  obvious  counterpart  in  die  estimated  Pn  and  Sn  anomalies,  although  higher  correlation 
is  observed  for  the  shear  wave  anomalies.  To  summarise  the  Norwegian  mountain  range 
correlates  poorly  with  crustal  thickness  and  does  not  appear  to  have  a  pronounced  root- 
structure  as  observed  elsewhere  in  similar  areas.  The  negative  P  and  S  velocity  anomalies 
observed  in  this  study  imply  mass  deficits  at  sub-Moho  levels  beneath  two  areas  of 
Norway.  In  this  context,  Meissner  (1979)  and  Theilen  and  Meissner  (1979)  have 
suggested  that  the  primary  mountain  roots  may  have  been  ’scraped  off  by  a  creep 
process  related  to  the  opening  of  the  Norwegian-Greenland  Sea.  The  implied  mechanism 
is  the  intrusion  of  low  viscosity,  low  density  asthenosphere  material  beneath  the 
Caledonides.  Such  a  process  is  likely  to  be  most  dominant  in  areas  close  to  the 
continental/  oceanic  transition,  i.e.,  close  to  where  the  velocity  anomalies  are  observed.  It 
is  also  notable  that  these  areas  (Mpre  and  Lofoten)  are  also  the  landward  extension  of 
the  major  fracture  zones  in  the  Norwegian  sea,  the  Jan  Mayen  Fracture  Zone  and  Senja 
Fracture  Zone  (Fig.9).  Realistically  however,  although  the  data  at  hand  indicates  the 
existence  of  sub-Moho  structural  heterogeneities,  the  data  are  not  sufficient  to  resolve  the 
type  and  extent  of  geodynamic  processes  deforming  the  lithosphere  in  areas  of 
W. Norway.  A  far  more  extensive  observational  data  base,  including  crustal  seismic 
structural  details  and  heatfiow  data,  is  required  for  a  serious  attempt  at  modelling  the 
collision  and  opening  processes. 
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The  geometry  of  station  and  event  distribution,  as  noted  above,  means  that  velocity 
estimates  within  pheripheral  blocks  are  less  well  resolved.  We  still  consider  however  that 
the  anomalies  coinciding  with  the  Fennoscandian  border  zone,  and  continuing  further 
westward,  reflect  sub-Moho  structural  heterogeneities  in  this  area.  The  anomalies  could  in 
fact  continue  southward  without  being  resolved  by  the  data  (see  Fig  6d).  The 

argument  here  is  that  the  crust  in  the  area  is  very  thin  (e.g.  see  Kinck  et  al,  1989)  and 
that  the  associated  relative  recent  deformation  processes  have  originated  in  the  deeper 
pan  of  the  lithosphere. 

7.  Conclusions 

Pn  and  Sn  velocity  estimates  are  derived  for  the  sub-Moho  beneath  Fennoscandia. 

Both  P-wave  and  shear  wave  velocity  anomalies  are  found  in  the  nonhwest  and 
southwest  of  Norway  and  in  central  Finland.  There  is  a  close  correlation  between  the 
anomalies  observed  in  Norway  and  observed  negative  Bouguer  gravity  anomalies, 
indicating  the  presence  of  spatially  restricted  regions  in  the  uppermost  mantle  involving 
low  (P  and  S)  velocity  and  corresponding  low  density.  Geologically  tlie  anomalies  are 
likely  to  be  associated  with  geodynamic  processes  tied  to  the  opening  of  the  Norwegian- 
Greenland  Sea. 

Acknowledgements 

This  research  was  supported  by  the  Royal  Norwegian  Council  for  Scientific  and  Industrial 
Research  (S.C.B.),  the  Defence  Advanced  Research  Projects  Agency  under  AFORS  Grant- 
89-0259  (B.O.R),  VISTA  (V6222)  and  NAVFCD  (E.S.H). 


77 


References 


Aki.K.  and  l^e,W.H.K.,1976.  Determination  of  three-dimensional  velocity  anomalies 
under  a  seismic  array  using  first  P  arrival  times  from  local  earthquakes  (part  1). 
J.geophys.Res.,81,  4381-4399. 

Aki,K.,  Christoffersson,A.  and  Husebye.E.S.,  1977.  Determination  of  the  three-dimensional 
seismic  structure  of  the  lithosphere.  J.Geophys.Res.,82;  277-296. 

Balling, N.,  1984.  Gravity  and  isostasy  in  the  Baltic  Shield.  In:  Galson,D.A.  and 
Mueller.S.  (eds),  First  EGT  Workshop  :  The  Northern  Segment,  53-68.  European  Science 
Foundation,  Strasbourg,  169  pp. 

Berthelsen,A.,  1987.  A  tectonic  model  for  the  crustal  evolution  of  the  Baltic  shield.  In  ; 
J.-P  Schaer  and  J.R.Rodgers  (eds.).  The  Anatomy  of  Mountain  Ranges.  Princeton 
University  Press.  Princeton,  New  Jersey,  31-57. 

Birch,F.,  1961.  The  velocity  of  compressional  waves  in  rocks  to  10  kilobars,  part  2. 
J.gcophys.Res,  66,  2199-2224. 

Calcagnile,G.  1982.  The  lithosphere-asthenosphere  system  in  Fennoscandia.  In: 

E.S.Husebye  (ed).  The  Structure  of  the  Lithosphere-Asthenosphere  in  Europe  and  the 
North  Atlantic.  Tectonophysics,  90:  19-35. 

Dcr,Z.A.  and  Landisman,M.1972.  Theory  of  errors,  resolution  and  separation  of  unknown 
variables  in  inverse  problems,  with  application  to  the  mantle  and  the  crust  in  southern 
Africa  and  Scandinavia.  Geophys.J.Roy.astron.Soc.,27,137-178. 

Dost,B.,  1986.  Preliminary  result  from  higher-mode  surface  wave  measurements  in 
western  Europe  using  NARS  array.  In:  D.A.Galson  and  St.Mueller  (eds.);  The  European 
Geotraverse,  Part  2.  Tectonophysics,  126,  1-30. 

Faleide,J.I.,  Gudlaugsson,S.T.,  Eldholin,0.,  Myhre,A.M.  and  Jackson,H.R.,  1989.  Deep  sea 
transects  across  the  western  Barents  Sea  continental  margin.  Tectonophysics,  this  issue. 

Gaal,  G.  and  Gorbatschev,R.  1987.  An  outline  of  the  Precambrian  evolution  of  the  Baltic 
shield.  Precambrian  Research,  35,  15-52. 

Galson,D.A.  &  Mueller,St.,  1986.  An  introduction  to  the  European  Geotraverse  project  : 
First  results  and  present  plans.  Tectonophysics,  126,  1-30. 

Gee,D.B.  and  Sturt, B. A.  (eds).  1985.  The  Caledonian  Orogen  -  Scandinavia  and  related 
areas.  Wiley,  New  York,  197-212. 

Mawonh,R.T.,  Hipkin,R.,  Jacobi,R.D.,  Kane,M.,  Lefort,J.P.,  Max,M.D.,  Miller,H.G.,  & 
F.Wolff.  1988.  Geophysical  framework  and  the  Appalachian-Caledonian  connection.  In: 
n:uTis,A.L.  &  Fcttes,D.j.  (eds),  7he  Caledonian-Appalachian  Orogen,  Geological  Society 
Special  Publication  No.  38,  pp  3-20. 

Hearn,  I. M.  and  Clayton, R.W.,  1986a.  Lateral  vekx'ity  variations  in  southern  California.  I. 
Results  for  the  upper  crust  from  Pg  waves.  Bull.Sei.sm.S(x;.Am.,76,  495-509. 


78 


Hearn, T-M.  and  Clayton, R.W.,  1986b.  Lateral  velocity  variations  in  southern  California.  2. 
Results  for  the  lower  crust  from  Pn  waves.  Bull. Seism. Soc.Am.,76,  511-520. 

HeiTnan,G.T.,  Lakshinmarayan,A.V.  and  Rowland;S.W.,  1975.  The  reconstruction  of 
objects  from  shadowgraphs  with  high  contrast.  Pattern  Rec.,  7;  157-165. 

Herrin,E.,1968.  Seismological  fables  for  P  phases.  Bull. Seism. Soc. Am.,  58,  1193-1241. 

Hestenes.M.R.  and  Stiefel,E.,1952.  Methods  of  conjugate  gradients  for  solving  linear 
systems.  J.Res.N.B.S., 49, 409-436. 

Hossack,J.R,  and  Cooper,M.A.,  1986.  Collision  tectonics  in  the  Scandinavian  Caledonides. 
In:  Coward,M.P.,  Dewey,J.F.  &  Hancock,P.L.(eds),  Continental  Extensional  Tectonics, 
Geological  Survey  Special  Publication  No.  28,  pp  611-619. 

Hovland,!.,  Gubbins,  D.,  and  Husebye,E.S.,  1981.  Upper  mantle  heterogeneities  beneath 
Central  Europe,  Geophys.J.R.astron.Soc.,  66,  261-284. 

Husebye.E.S.,  1982.  (Editor),  The  structure  of  the  lithosphere-asthenosphere  in  Europe  and 
the  North  Atlantic.  Tectonophysics,  90(l-2):l-262. 

Husebye,E.S.  and  HovlandJ.,1982.  On  upper  mantle  seismic  heterogeneities  beneath 
Fennoscandia.  Tectonophysics,  90:  1-17. 

Husebye,E.S.,  Hovland,!.,  Christoffersson,A.,  Astrom,K.,  Slunga,R.  and  C-E.Lund,  1986. 
Tomographic  mapping  of  the  lithosphere  and  asthenosphere  beneath  southern  scandanavia 
and  adjacent  areas.  Tectonophysics,  128,  229-250. 

Kinck,J.J.,  Husebye,E.S.  and  C.-E.  Lund.,  1989.  The  S.  Scandinavia  Crust  -  Structural 
complexities  from  seismic  reflection  and  refraction  profiles.  Tectonophysics,  this  issue. 

KnopoffL..  1972.  Observation  and  inversion  of  surface  wave  dispersion.  Tectonophysics 
13:  497-513. 

Luosto.U.  1989.  Crustal  structures  of  eastern  Fennoscandia.  Tectonophysics,  this  issue. 

Meissner,R.  1979.  Fennoscandia  -  a  short  outline  of  its  geodynamical  development.  Geo 
Journal,  3.3,  227-233. 

Nolet,G.,  1985.  Solving  or  resolving  inadequate  and  noisy  tomographic  systems, 
J.Comp.Phys.,  61,  463-482. 

Paige,C.C.  and  Saunders, M. A. ,1982.  LSQR:  an  algorithm  for  sparse  linear  equations  and 
sparse  least  squares,  ACM  Trans. Math.Soft.,8,  43-71  and  195-209. 

Planke,S.,  Skogseid,J.  and  Eldhom,0.,  1989.  Crustal  structure  off  Norway,  62°  to  70° 
North.  Tectonophysics,  this  issue. 

ProdchI,C.  and  Kaminski, W.,  1984.  Crustal  structure  under  the  FENNOLORA  profile.  In  : 
Galson,D.A.  and  MueIler,S.  (eds).  First  EGT  Workshop  'i  ne  Northern  Segment,  53-68. 
European  Science  Foundation,  Strasbourg,  169  pp. 

Ruud.B.O.,  Husebye,E.S.,  Ingate,S.F.,  and  Chrisvuk  r.sS'ai,A.,  1988.  Event  location  at  any 


79 


distance  using  seismic  data  from  a  single,  three-component  station,  Bull.Seism.Soc.Am. 
78,  308-315. 

Spakman,W.  and  G.Nolet,  1988.  Imaging  algorithms,  accuracy  and  resolution  in  delay 
time  tomography,  in.  Mathematical  Geophysics;  a  survey  of  recent  developments  in 
seismology  and  geodynamics,  edited  by  N.J.Vlaar  et  al.,  Reidel,  Dordrecht,  155-188. 

Spencer,C.P.  and  Gubbins,D.  1980.  Travel-time  inversion  for  simultaneous  earthquake 
location  and  velocity  structure  determination  in  laterally  varying  media, 
Geophys.J.R.astr.Soc.,63,  95-116. 

Talwani,M.  and  Eldholm,0.  1977.  Evolution  of  the  Norwegian-Greenland  Sea. 
Geol.Soc.Am.Bull.  88,  969-999. 

Theilen,F.  and  Meissner,R.  1979.  A  comparison  of  crustal  and  upper  mantle  features  in 
Fennoscandia  and  the  Rhenish  Shield,  two  areas  of  recent  uplift.  Tectonophysics,  61,  227- 
242. 

Wolff,F.  1983.  Crustal  structure  of  the  Scandinavian  Peninsula  as  deduced  by 
wavelength  filtering  of  gravity  data.  In:  Schenk,P.E.  (ed.)  Regional  Trends  in  the 
geology  of  the  Appalachian-Caledonian-Hercynian-Mauritanide  Orogen,  Reodel, 

Dordrecht. 

Ziegler,P.A.  1988.  Evolution  of  the  Arctic-North  Atlantic  and  the  Western  Tethys. 

AAPG  Memoir  44,  USA. 


80 


Figure  captions; 

Fig.l;  The  model  of  Moho  depth.  The  shading  is  from  9  km  (dark)  to  55  km  thickness 
(light  shade). 

Fig.2:  The  distribution  of  events  (circles)  and  seismic  stations  (triangles)  within 
Fenncscandia.  Circled  symbols  refer  to  the  national  symbols  of  Norway  (N),  Sweden 
(S),  Finland  (SF)  and  Denmark  (DK).  P.Z.  =  protogine  zone,  B.B.  =  Bothnian  Bay, 
S.F.Z.  =  Senja  Fracture  Zone,  JMFZ  =  Jan  Mayen  Fracture  Zone. 

Fig.3:  Reduced  travel  time  (V,^  =  9  km/s)  versus  distance  of  all  arrivals  reported  in  the 
Fennoscandian  seismic  data  base,  Jan  1984  -  May  1988.  The  polygon  windows  delineate 
the  Pn  (lower  window)  and  Sn  (upper)  arrivals  that  were  accepted  in  the  initial  data 
selection. 

Fig.4:  The  spatial  coverage  of  Pn  raypaths  used  in  the  inversion.  The  Sn  raypaths  is  not 
visually  different. 

Fig.5:  Epicenters  of  the  events.  The  circles  refer  to  the  epicentral  locations  following 
the  inversion  procedure,  the  attached  lines  indicating  the  shift  from  the  initial  Finland 
data  center  locations. 

Fig.6:  Reconstructions  from  synthetic  data,  (a)  Synthetic  model  involving  a  pattern 
anomaly  with  blocks  of  velocity  7.9  km/s  (dark  shade)  and  8.5  km/s  (light  shade)  set 
within  a  background  velocity  of  8.2  km/s.  (b)  The  reconstruction  of  this  synthetic 
model,  shading  ranges  from  8  km/s  (dark)  to  8.36  km/s  (light),  (c)  Synthetic  model 
involving  horizontally  oriented  anomalies,  initial  velocities  as  in  (a),  (d)  Reconstruction 
of  (c),  velocity  range:  8  km/s  (dark)  -  8.3  km/s  (light). 

Rg.7:  P-wave  velocities  estimated  from  the  real  data.  Shading  ranges  from  8.04  km/s 
'ark)  to  8.3  km/s  (light).  A  reference  velocity  of  8.18  km/s  was  used  in  the  initial 
i  lei. 

Fig. 8:  S-wave  velocity  estimates.  Shading  from  4.61  km/s  (dark)  to  4.69  km/s  (light).  A 
reference  velocity  of  4.65  was  used  in  the  initial  model. 

Fig.9:  Geology  of  Fennoscandia  (following  Berthelsen,  1987).  SF  -  Svecofennian 
province,  PZ  -  Protogine  Zone,  SN  -  Sveconorwegian  orogen,  A  -  Archean  nucleus, 
JMFZ  -  Jan  Mayen  Fracture  Zone,  SFZ  -  Senja  Fracture  Zone,  ESC  -  Faeroe- Shetland 
and  Vdring  Plateau  escarpments  (southern  and  northern  ESC  respectively),  FBZ  - 
Fennoscandian  Border  Zone.  Bathymetry  is  marked  in  hundreds  of  meters  (dotted  lines). 

Fig.  10:  Map  of  long  wavelength  Bouguer  gravity  anomalies  (Balling,  1984).  Contours 
are  at  20  mgal  intervals. 
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Abstract 

P-wave  coda  of  75  s  duration  from  8  teleseismic  events  of  widely  different  azimuths 
and  depths  as  recorded  by  the  NORESS  array,  Norway,  have  been  examined  using  both 
array  and  3C  analysis  techniques.  The  array  data  have  been  analyzed  using  the 
semblance  technique  which  allow  us  to  determine  with  good  accuracy  both  time  of 
arrival,  apparent  velocity  and  azimuth  of  the  scattered  waves.  Coda  coherency  was 
found  to  vary  considerably,  probably  reflecting  extended  source  duration  for  some 
events.  By  subtracting  the  teleseismic  beam  from  the  records,  locally  scattered  wavelets 
can  be  resolved  even  within  highly  coherent  teleseismic  P-coda.  Although  source-end 
scattering  could  not  be  separated  from  the  source  pulse  per  se,  we  favor  the  hypothesis 
of  long  source  duration  supported  by  observed  slow  beam  amplitude  decay  rates.  The 
majority  of  receiver-end  scattering  contributions  appears  to  be  P-to-Rg  conversion  in 
both  forward  and  backward  modes  from  two  nearby  areas  with  pronounced  topographic 
reliefs,  namely  Bronkeberget  (Dist~10  km,  Azi~80°)  and  Skreikampen  (Dist~30  km, 
Azi~225°).  The  scattering  is  multiple  in  the  sense  that  both  the  primary  and  secondary 
phases  from  the  source  region  contribute  to  rhe  Rg-scattering.  P  to-S  scattering 
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constitute  a  significant  part  of  receiver-end  scattering  but  is  more  diffuse  than  Rg.  Also 
here  the  Rg  preferred  azimuth  directions  am  found,  but  the  forward  direction  is 
important  as  well.  P-to-P  scattering  is  weak  and  mostly  confined  to  the  immediate 
vicinity  of  the  array.  The  analysis  of  individual  3C  stations  shows  high  sensibility  to 
interference  with  locally  scattered  waves  which  are  probably  responsible  for  a  marked 
decrease  in  apparent  velocity  observed  for  some  events  3-4  s  after  P-onset.  To  fully 
utilize  the  potential  of  3C  instruments  for  phase  identification  in  an  integrated  array/3C 
processing  scheme,  at  least  one  additional  3C  station  within  the  NORESS  array  would 
be  necessary. 


1.  Introduction 

The  classical  earth  model  involving  spherically  symmetric  layers  is  valid  only  as  a  first 
approximation  besides  being  useful  in  the  computation  of  synthetic  seismograms. 
Although  dominant  and  deterministic  seismogram  features  like  P,  S  and  surface  waves 
can  be  accounted  for  using  such  models,  their  codas  cannot  be.  Coda  waves  arriving 
after  deterministic  phases  are  broadly  attributed  to  seismic  wave  scattering  phenomena 
associated  with  lateral  heterogeneities,  located  within  the  crust  and  mantle  while 
precursors  to  PKIKP  and  related  phases  reflect  heterogeneities  near  the  core-mantle 
boundary.  Velocity  and  density  heterogeneities,  in  addition  to  generating  coda  waves, 
can  cause  changes  in  waveform,  phase  (or  travel  time)  and  amplitude  fluctuations  as 
well  as  attenuation  of  the  direct  (deterministic)  phas^'  arrivals.  Overall  the  seismic  wave 
scattering  reflects  the  dominant  seismic  wavelength  of  he  seismic  signal  relative  to  the 
scale  length  of  the  heterogeneities,  which  may  vary  from  the  grain  size  of  rocks  to  the 
lowest  mode  of  global  spherical  harmonics;  an  instructive  overview  is  given  by  Wu  and 
Aki  (1988) 

Seismic  wave  scattering  in  a  heterogeneous  earth  has  motivated  numerous  theoretical 
studies.  Most  of  these  studies  are  aimed  at  obtaining  approximate  solutions  to  the  wave 
equation  for  a  medium  characterized  by  a  few  parameters  such  as  scale  length  and 
perturbation  index  (rms  variation  of  velocity  and  density).  Other  studies  concentrate  on 
calculating  die  seismic  response  of  wedge  suuctuics,  Moho  and  topographic  undulations 
and  geological  structures  of  irregular  geometries.  In  this  kind  of  studies  it  is  not  always 
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clear  which  propagation  effects  (e.g.  mode  conversion,  diffraction  and  medium-related 
anisotropy)  are  included.  Some  of  these  problems  can  be  reduced  by  using  finite 
element  methods  (e.g.,  see  Frankel,  1989  and  Herraiz  and  Espinosa,  1987). 

Observational  studies  of  wave  scattering  effects  are  problematic  due  the  weak  signal 
correlation,  even  when  using  small  arrays.  This  seems  to  be  one  reason  why 
observational  scattering  parameters  are  often  tied  to  coda  envelopes  and  decay  rates 
while  wave  types  and  associated  propagation  paths  constituting  the  coda  waves  are  often 
not  specified.  An  exception  is  naturally  the  amplitude  and  phase  fluctuations  of  initial 
P-wave  arrivals,  as  demonstrated  by  Aki  (1973)  and  Berteussen  et  al  (1975)  using  the 
random  media  theory  of  Chernov  (1960).  These  studies  led  to  the  tomographic  mapping 
approach  or  the  ACH-method  of  Aid  et  al  (1977). 

In  this  study  we  focus  on  attributes  of  P-coda  waves  as  recorded  by  the  small  aperture 
NORESS  array  in  southern  Norway  (Fig.  1  and  Ingate  et  al,  1985).  The  aims  of  our 
study  are  two-fold,  namely  to  attempt  to  decompose  the  coda  waves  into  specific  phase 
arrivals  and  associated  slowness  vectors  and  secondly  to  provide  an  insight  into  the 
scattering  processes  most  important  in  generating  the  coda  waves  of  teleseismic  P- 
signals. 


2,  P  coda  generation  -  source  or  receiver  contributions  ? 

An  important  parameter  in  scattering  theory  is  the  dimensionless  quantity  ka  where  k  is 
the  wavenumber  and  a  the  correlation  distance  in  random  media  theory  or  scale  length 
of  the  structural  heterogeneities.  For  teleseismic  P-coda  with  frequencies  in  the  range  1- 
3  Hz  and  structural  heterogeneities  of  the  order  of  0.1  to  10  km,  the  factor  ka  may 
vary  from  0.1  to  ca  30.  The  latter  sizes  of  heterogeneities  have  been  confirmed  from  P- 
phase  and  amplitude  studies  in  the  area  of  die  NORESS  site  (e.g.  see  Aki  et  al.,  1977, 
Haddon  and  Husebye,  1979  and  Flatt6  and  Wu,  1988)  while  the  smaller  scale  length  of 
ca.  1  km  or  less  has  been  observed  in  many  scattering  studies  around  the  world  from 
coda  analysis  of  high  frequency  local  events.  Small  angle  or  forward  scattering  would 
be  important  for  ka  .>.>  1  in  which  case  sections  of  the  coda  should  realistically  be  in 
the  source  direction.  When  ka  ~  1  however,  wide  angle  scattering  contributions  are 
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expected  to  be  relatively  strong,  with  large  angle  scattering  relative  to  a  given  azimuth 
direction.  Conversions  between  different  wave  types  are  by  far  most  efficient  when 
structural  anomalies  are  both  in  velocity  and  density.  The  crust  and  the  surface 
topography  are  on  this  basis  likely  to  play  an  important  role  in  scattering  generation. 
Second-order  or  doubly  scattered  wavelets  are  generally  considered  to  be  of  little 
importance  since  the  earth  appears  to  be  weakly  heterogeneous,  i.e.  density  and  velocity 
rms  fluctuations  seldom  exceed  5%.  The  relative  excitation  of  the  coda  may  vary  with 
wavelength  or  equivalent  signal  frequencies  given  that  the  scatter  length  distribution  is 
peaked,  say,  around  10  km. 

Scattering  may  conceptionally  be  thought  of  as  occurring  in  three  principal  areas, 
namely  in  the  source  region  of  the  earthquake,  along  the  wave  path  in  the  mantle  (say 
below  4(X)  km)  and  finally,  near  the  receiver.  We  discuss  each  of  these  separately 
below. 

Source  area  scattering  :  The  size  and  focal  depth  of  an  eanhquake  is  obviously  of 
some  importance,  at  least  in  regards  to  the  time  duration  of  scattering  within  the  source 
region.  Scattering  contributions  may  arise  from  either  multiple  P-reflections  within  the 
crustal  and  lithospheric  layers  and/or  S-  or  R-wave  conversions  to  P.  Intuitively  the 
coda  contribution  of  near-source  scattering  is  expected  to  be  relatively  coherent  when 
observed  at  teleseismic  distances  because  of  the  small  variation  in  take-off  angle.  To 
study  source-end  scattering  it  would  be  necessary  to  use  data  from  several  stations,  e.g. 
see  Lynnes  and  Lay  (1989). 

Transmission  path  scattering  :  For  relatively  short  distances  of  20-30°  layered 
structures  within  the  lithosphere  and  asthenospherc  may  act  as  a  waveguide,  thus 
contributing  significantly  to  the  P  coda  as  demonstrated  by  Kennett  (1987)  and  Kom 
(1988).  Somewhat  similar  effects  have  been  noted  by  King  et  al  (1977)  in  the  analysis 
of  PP  precursors  where  distinct  differences  between  observed  slownesses  were  attributed 
to  P400PP  at  the  source-end  and  PP400P  at  the  receiver  end  of  the  propagation  path. 

Receiver  end  scattering  :  In  general  there  is  some  differentiation  here  between  two 
effects;  that  resulting  from  large  scale  heterogeneities,  including  continental  and  oceanic 
transition  structural  regimes,  such  as  off  the  west  coast  of  Norway,  and  secondly,  small 
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scale  scatterers  in  the  crust  and  lithosphere,  including  rough  topography.  P-to-Rg 
conversion  should  be  effective  in  the  latter  case  if  the  topography  undulations  at  least 
are  of  the  order  of  half  the  Rayleigh  wavelength  of  approximately  1-2  km  (Frankel, 
1989).  The  west  Norway  mountain  area  ca  100-300  km  away  from  the  NORESS 
seismic  array  is  thus  a  candidate  for  this  kind  of  scattering. 


3.  Data  and  analysis  technique 

As  several  coda  wave  mechanism  may  contribute  we  attempt  to  isolate  individual 
scattering  sources  through  careful  event  selection.  Table  1  indicates  good  azimuth 
coverage,  thus  illuminating  potential  scatterers  from  a  multitude  of  directions.  We  have 
used  both  3C  recordings  from  the  four  3C  stations  within  the  NORESS  array  (Fig.  1) 
and  vertical  component  data  from  the  whole  array.  The  essential  signal  parameters 
extracted  for  deducing  wave  type  and  propagation  paths  are  polarity  (particle  motion) 
and  slowness  vectors  in  specific  time  windows. 

Analysis  techniques 

In  coda  analysis  optimum  use  of  the  various  filtering  options  in  the  context  of  array 
wave  field  recording  is  important.  To  suppress  pan  of  the  ambient  noise  band-bass 
filtering  in  the  range  0.5-3.0  Hz  are  used.  Source-end  scattering  contributions  and  the 
signal  pulse  itself  are  taken  to  be  equivalent  to  the  beam  trace  resulting  from  aiming 
the  array  towards  the  epicenter.  The  residual  trace,  that  is,  the  difference  between  single 
sensor  records  and  the  beam  trace,  would  enhance  receiver-end  scattering  contributions. 
Beam-forming  on  the  basis  of  residual  traces  would  then  give  a  relative  powerful  tool 
for  pin-pointing  scattered  wavelets.  Simple  rotation  of  the  NS-  and  EW-components  into 
radial  and  transverse  components  relative  to  the  epicenter  and/or  scattering  sources 
would  naturally  serve  to  isolate  specific  wave  types  in  the  records. 

Data  analysis  included  f-k  and  semblance  analysis  (Husebye  and  Ruud,  1989),  both  of 
which  have  good  separation  of  body  and  surface  waves  even  for  a  small  aperture  array 
like  NORESS.  We  also  experimented  with  semblance  velocity-analysis  of  the  radial  and 
transverse  components  of  the  four  3C  stations  (Fig.  1),  but  the  results  were  not 
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instructive  due  to  the  odd  configuration  of  Lhis  mini-array.  At  least  one  additional  3C 
station  is  needed  in  the  C-ring. 

3C  analysis  techniques  available  include  that  of  Vidale(1986),  Park  et  al  (1987),  Plesinger 
et  al  (1986),  Jurkevics  (1988),  Christoffersson  et  al  (1988)  and  Roberts  and 
Chi  stoffersson  (1990).  We  have  used  the  latter  in  our  coda  analysis.  This  technique 
works  in  the  complex  demodulate  domain,  a  complex  time  domain,  efficiently 
decomposing  3C  signals  to  obtain  instantaneous  estimates  of  the  real  and  quadrature  parts 
(r  nplitude  and  phase)  of  the  energy  within  selected  frequency  bands.  Various  attributes 
of  the  seismic  coda  can  be  easily  obtained  from  the  array  and  3C  records.  Below  we 
outline  several  attributes  used  in  this  study. 

Wavelet  azimuth;  The  azimuth  is  estimated  from  the  observed  cross-spectral  densities 
between  the  horizontal  and  vertical  components.  Although  this  estimate  is  instantaneous, 
some  averaging  within  a  window  of  around  1  cycle  is  usually  done  to  improve  stability. 

In  the  analysis  of  array  data  the  slowness  vector  (apparent  velocity  and  azimuth)  is  tied 
to  the  maximum  beam  power  (f-k)  or  the  maximum  coherence  (semblance). 

Incidence  angle  and  apparent  velocity:  The  incidence  angle  is  determined  from  the 
cross-spectral  density  between  the  vertical  and  radial  components.  This  angle  can  be 
converted  to  a  velocity  estimate  given  assumptions  about  the  local  P-wave  velocity  and 
Poissons  ratio  or  using  known  epicenter  locations  for  calibration.  For  an  array  the 
apparent  velocity  is  estimated  directly  from  move-out  times;  the  resolution  is  poor  when 
the  apparent  wavelength  of  the  signal  greatly  exceeds  the  array  aperture,  and  for 
NORESS  this  apply  to  the  entire  teleseismic  window. 

Phase  difference  between  components:  We  can  examine  the  instantaneous  phase 
difference  between  the  rotated  radial  and  vertical  component,  or  equally  between  other 
pairs  of  components,  using  the  phase  information  retained  in  the  3-component  analysis.  A 
phase  difference  of  180°  corresponds  to  pure  P  wave  particle  motion,  0°  difference 
represents  pure  SV  wave  particle  motion  while  a  phase  difference  of  ±90°  represents 
elliptical  particle  motion,  as  expected  for  example  for  Rayleigh  waves  but  also  for  SV 
waves  with  post-critical  incident  angles. 
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Array  semblance  and  3C  signal  polarity:  Our  confidence  in  the  above  derived  signal 
attributes  depends  on  the  associated  semblance  and  polarity  estimates.  Semblance  values 
(AO  plus  C-  and  D-ring)  for  random  data  seldom  exceeds  0.2  for  the  frequency  range  and 
time  window  lengths  used.  If  coi responding  values  are  observed  in  the  coda  sections  they 
are  considered  marginal  or  not  significant  unless  the  associated  slowness  vector  estimates 
remain  stationary  over  a  few  cycles.  For  3C  analysis  we  have  different  measures  on  the 
extent  of  polarity  of  the  signals  at  hand,  including  the  chi-square  model  fit  measure  of 
Christoffersson  et  al  (1988)  where  the  significance  is  marginal  for  values  below  0.5.  To 
avoid  lengthy  explanations  we  use  a  common  denominator  for  both  semblance  and 
polarity  measures  -  the  term  signal  coherency. 

Practical  considerations:  The  major  difficulty  faced  in  coda  analysis  is  that  of  dealing 
with  interference  resulting  from  multiple  wavelets  arriving,  at  least  partly,  simultaneously. 
In  such  a  case  the  recorded  signal  involves  rapid  variations  in  the  polarization,  with 
fluctuations  in  both  calculated  azimuth  and  angle  of  incidence.The  above  signal  attributes 
are  calculated  using  a  sliding-window  approach.  Window  lengths  are  around  1-2  cycles 
(1-2  s)  while  the  updating  frequency  is  1/3  or  1/4  of  this  length.  Overall  we  concentrate 
on  separating  signals  where  the  polarization,  phase,  semblance  and  slowness  are  such  that 
clearly  definable  phases  exist,  examining  their  consistency  over  time  and  their  spatial 
correlation  across  the  NORESS  array. 


4.  Results 

Altogether  8  events  were  subject  to  analysis  and  the  results  would  be  presented  along  the 
classical  scattering  subdivision  of  source-  and  receiver-end  contributions.  Then  appropriate 
individual  event  results  would  be  addressed  in  some  detail,  Hypocenter  information  is 
listed  in  Table  1  while  trace  displays  are  in  Fig.  2.  In  general  we  have  that  events  to  the 
south  and  west,  that  is  MAD,  LEE,  ALA  (Fig.  2b)  are  of  relatively  low  frequency  and 
besides  these  events  also  have  a  relative  long  P-signal  duration.  The  latter  also  applies  to 
the  HKl  event  (Fig.  2a),  a  Hindu  Kush  earthquake  with  a  focal  depth  of  100  km,  which 
appear  to  be  a  double/triple  event.  A  quantitative  measure  of  P-signal  duration  including 
source-end  scattering  is  given  in  Table  2  and  Fig.  3.  Note  that  the  deep  events  JAP  and 
HK2  and  the  SEM  explosion  have  a  very  strong  source  pulse  of  3-5  s  duration  (Fig  2a), 
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but  the  source-end  scatteri  ig  contributions  result  in  a  total  duration  of  more  than  60  s  in 
terms  of  coherent  P-telcseismic  arrivals  (Fig.  3). 

4.1  Source-end  scattering. 

I'he  prime  characteristics  here  would  be  that  scattered  pulses  are  P-type  when  observed  at 
the  receiver  in  the  early  P-coda  and  this  naturally  implies  that  scattering  and  complex 
source  pulses  arc  not  easily  separated.  Anyway,  3C  techniques  were  found  convenient  for 
source-end  scattering  analysis,  and  relevant  results  here  are  displayed  in  Fig.  2.  To 
familiarize  the  reader  with  the  3C-results  these  are  discussed  in  some  detail  for  the 
Kamchatka  (KAM)  event  (Fig.  2a).  The  records  for  site  AO  have  been  bandpass  filtered 
and  rotated  with  the  radial  component  in  the  source  azimuth  direction.  Processing 
parameters  are  explained  in  the  figure  caption.  An  important  note  is  that  the  azimuth- 
coherency  and  velocity  plots  in  the  Fig.  b  and  d  are  based  on  the  assumption  of  P-wave 
presence  while  there  is  no  such  assumption  on  wave-type  with  the  calculated  phase  shift 
between  vertical  and  radial  components.  Note  that  for  wavelets  with  azimuth  directions 
greatly  off  ’true’  azimuth  the  observed  phase  shifts  are  not  directly  interpretable  in  P,  S 
and  Rg-types.  Signal  coherency  and  hence  confidence  in  azimuth  estimates  are  good  for 
the  first  7  s  of  the  signal  (5-12  s)  as  seen  in  the  upper  left  corner  of  the  figure.  A  small 
discontinuity  at  8-9  s  has  a  counterpart  in  the  R-Z  phase  diagram  and  besides  marks  a 
transition  to  a  pronounced  drop  in  apparent  velocity.  Between  12-13.5  s  signal  coherency 
is  low  (below  0.5)  while  the  phase  shift  drops  toward  90°,  implying  elliptical  particle 
wave  motion.  From  13.5  to  15  s  the  signal  coherency  is  again  high  and  the  associated 
phase  shifts  are  typical  P-wave.  This  contrasts  with  the  values  between  17  and  20  s 
where  coherency  is  good  but  the  phase  shifts  in  Fig.  3c  imply  elliptical  motion.  In  the 
remaining  parts  of  the  displayed  records,  signal  coherency  is  mostly  low  and  particle 
motion  is  predominantly  elliptical,  fluctuating  mainly  between  90  and  270°. 

Corresponding  velocity  estimates  are  calculated  at  a  lower  coherency  level  (above  0.2), 
subject  to  the  condition  that  the  standard  deviation  of  the  estimate  is  less  than  1.6  of  the 
actual  velocity  value.  The  very  reliability  of  the  obtained  velocity  estimate  is  probably 
poor,  the  low  coherency  representing  strong  interference  from  non-P  wavelets.  To 
summarize  the  3C  results  for  the  KAM  event,  strongly  coherent  P-wave  motion  in  the 
first  part  of  the  record  together  with  the  estimated  azimuth  and  apparent  velocity  are  very 
close  to  that  expected  for  the  event  location.  After  about  3  s  we  have  a  sharp  drop  in 
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velocity  although  the  azimuth  does  not  change  significantly. 

The  3C  results  from  the  other  7  events  also  in  Fig.  2  are  not  much  different  from  those 
typical  for  KAM.  The  common  characteristics  are  as  follows;  the  phase  ihift  results 
imply  that  interference  starts  3-4  s  after  P-onset  and  often  are  so  prominent  that  velocity 
and  azimuth  estimation  fails.  Also,  phase  shifts  fluctuate  between  90-270°  and  thus 
indicate  a  dominance  of  elliptical  (Rayleigh)  wave  motion.  Later  arriving  P-wavelets  are 
seldom  although  the  HKl,  MAD,  LEE  and  ALA  events  have  extended  P-signal  duration. 
Another  feature  is  that  the  P-velocity  sometimes  drop  significantly  after  just  3-5  s  (e.g. 
KAM,  LEE,  ALA).  This  is  rather  puzzling  since  the  time  lag  between  direct  P  and 
’secondary’  P-wavelets  are  small,  and  thus  difficult  to  explain  except  for  sort  of 
multipathing  confined  to  the  azimuth  plane.  Using  the  array  per  se,  we  have  calculated 
apparent  velocities  using  AO,  C-,  and  D-ring  and  here  the  observed  velocity  decreases  are 
very  moderate.  However,  the  3C  results  have  a  counterpart  in  the  displayed  traces  as  the 
energy  on  the  radial  components  increase  faster  than  the  vertical  components  after  the  P- 
onset.  From  Table  2  it  is  obvious  that  the  vertical  components  are  rather  coherent  across 
the  array  compared  to  the  3C  records  being  far  more  sensitive  to  interference.  This  point 
is  further  illustrated  in  Fig.  3. 


4.2  Receiver-end  scattering. 

From  the  above  discussion  and  also  from  Fig.  3  it  is  clear  that  source-end  contributions 
are  rather  energetic.  We  have  removed  such  effects  by  forming  the  beam  trace, 
subtracting  this  from  the  single  z-traces  and  then  conducting  semblance  analysis  on  the 
residual  traces.  We  have  screened  the  coda  in  two  segments  of  40  s  (Fig.  2  shows  the 
traces  in  the  first  segment)  for  P  (6-8  km/s),  Sn  (4.5  km/s)  and  Rg  (3.0  km/s).  In  the 
first  time  segment,  the  coda  appears  to  be  dominated  by  Rg-contributions  as  shown  i  Fig. 
4.  The  truly  striking  features  are  that  the  P-to-Rg  scattered  wavelets  mainly  come  from 
two  principal  azimuth  directions  of  80°  and  225°  (  thin  lines  in  the  figure)  irrespective  of 
event  azimuths.  Examples  of  estimated  Rg  waveforms  are  shown  in  Fig.  5.  With  the  time 
lags  involved,  we  have  located  the  corresponding  scattering  source  areas  to  Bronkeberget 
(Dist~10  km,  Azi~80°)  and  Skreikampen  (Dist~30  km,  Azi~225°).  These  areas  have 
prominent  topographic  reliefs,  and  are  marked  in  Fig.  1,  Furthermore,  the  scattering 
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appears  to  be  multiple  in  the  sense  that  Rj-wavelets  arrive  repeatedly  from  the  two 
azimuth  directions.  This  is  perhaps  most  clearly  illustrated  for  the  ALA  event  where  the 
timelags  between  scattering  pulses  on  the  80°  and  225°  azimuths  forms  a  distinct  pattern. 
Seemingly  Rg-attenuation  is  very  strong  since  wavelets  from  larger  distances  and  other 
azimuths  are  hardly  observed,  not  at  least  in  view  of  mountainous  area  to  the  north  and 
west  (60  km  or  more  away).  Semblance  analysis  was  performed  for  4.5  km/s  velocity 
and  the  main  results  here  (not  shown)  have  many  similarities  to  those  in  Fig.  4.  With 
this  is  meant  that  the  two  semblance  patterns  have  a  considerable  overlap  probably  due  to 
true’  velocities  around  3. 6-3.8  km/s  typical  of  Sg-phases.  However,  from  the  phase  shift 
results  and  from  a  good  overlap  in  azimuth  and  time  we  consider  Rg-wavelets  to 
dominate.  An  exception  here  is  the  high  semblances  for  ALA  at  about  50°  (30-35  s) 
which  are  much  stronger  in  the  4.5  km/s  results.  A  search  for  crustal  P-phases  (6-8  krr/s 
velocity  range)  turned  out  to  be  mainly  negative.  We  found  a  few  Pg-phases  at  4-6  s 
after  P-onset  time  but  hardly  later  in  the  coda.  This  implies  that  P-to-P  scattering  in  the 
general  array  siting  area  is  weak. 

Results  from  semblance  analysis  in  the  second  time  segment  (35-75  s  after  P-onset)  for 
velocities  of  3.0  to  4.5  km/s  are  shown  in  Figs.  6  and  7.  There  are  still  some  Rg- 
contributions  apparent  and  again  mainly  along  the  80°  and  225°  azimuth  directions.  The 
preference  for  the  above  two  azimuth  directions  in  Fig.  6  we  again  take  to  signify  that 
the  scattering  is  multiple  since  source -end  arrivals  are  still  present  far  into  the  coda  as 
shown  in  Fig.  3.  From  ’'ig.  7  we  see  tliat  S-type  scattering  dominates  in  this  time 
segment.  Tnis  is  explained  in  terms  of  relative  less  pronounced  S-wave  attenuation  as 
compared  to  Rg,  as  over  time  more  distant  scattering  source  areas  may  become  activated. 
Still  there  is  some  preference  for  the  80°  and  225°  azimuth  directions,  but  also  forward 
scattering  appears  to  be  significant,  in  particular  for  the  two  northern  events  KAM  and 
ALA,  but  also  for  SEM  and  HK  events  to  the  east.  P-to-P  scattering  at  the  receiver-end, 
tied  to  crustal  velocities  in  the  range  6-8  km/s,  remains  weak. 

5.  Discussion  and  Conclusions 

The  P-teleseismic  coda  has  also  been  examines  by  others,  and  we  would  briefly  review 
the  outcome  of  such  studies.  Already  two  decades  ago  Key  (1967)  identified  discrete 
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arrivals  in  the  P-coda  at  the  Eskdalemuir  array  in  Scotland  wliicn  suggested  P  to 
Rayleigh  wave  scattering  stemming  from  the  local  topography.  Recently  Langston  (1989) 
inferred  Rayleigh  wave  scattering  just  after  the  P-arrival  in  3C  analysis  of  teleseismic 
recordings  at  Pasadena,  California.  Of  most  interest  for  our  study  is  invesngations  by 
Dainty  and  Harris  (1989),  Dainty  (1990)  and  Gupta  et  al  (1989)  on  the  basis  of  f-k 
analysis  of  NORESS  teleseismic  codas.  Dainty  and  Harris  found  evidence  of  significant 
low-velocity  scattering  contributions  mostly  attributed  to  Lg  and  surface  waves  without 
detailing  azimuth  or  time  dependencies.  Gupta  et  al  reported  low  velocity  secondary 
arrivals  in  the  early  P-coda  to  the  east  and  south-west.  Our  preference  for  using  the 
semblance  technique  in  the  array  data  analysis  instead  of  high  resolution  f-k  methods  is 
due  to  the  fact  that  semblance  works  well  even  for  short  windows  (1-2  cycles)  and  thus 
the  individual  slowness  estimates  can  be  directly  related  to  arrivals  seen  in  the  records. 

As  also  demonstrated  here,  the  3C  records  are  useful  for  pin-pointing  where  wave 
interference  take  place,  and  naturally  the  information  potential  of  the  horizontal 
components  should  be  realized  in  coda  studies. 

5,1  Mechanism  for  NORESS  P-coda  generation: 

Although  teleseismic  scattering  is  generally  rated  a  complex  problem,  the  NORESS  coda 
observations  as  presented  in  the  previous  section  can  be  explained  fairly  simply.  Firstly, 
we  consider  source-end  contributions  to  be  of  long  duration  (more  than  one  minute)  even 
for  the  deep  events  like  JAP,  HK2  and  the  SEM.  At  the  source  end  a  number  of 
wavelets  are  generated,  the  dominant  ones  stemming  from  the  source  per  se  while  the 
others  are  secondary  ones  related  to  reverberations,  reflections,  mode  conversions  etc  in 
the  source  region  (see  also  Lynnes  and  Lay,  1989).  The  common  denominator  here  is 
apparent  velocities  appropriate  for  teleseismic  ranges.  At  the  receiver-end,  P-to-Rg 
conversion  related  to  topographic  relief  appears  to  be  most  efficient.  Naturally,  the 
scattering  areas  cannot  be  too  far  away  due  to  strong  Rg  attenuation  -  in  case  of 
NORESS  Rg  is  seldom  observed  beyond  100  km  for  local  events.  Anyway,  roughly  the 
first  30  s  of  the  coda  is  dominated  by  Rg-wavelets  originating  from  the  Bronkeberget  and 
Skreikampen  areas.  The  scattering  ’illumination’  of  NORESS  appear  to  be  dependent  on 
event  azimuth  and  is  quite  weak  for  some  events  like  LEE  and  ALA.  Although  the  two 
mentioned  azimuth  directions  dominates  for  most  of  the  events,  the  efficiency  of 
particular  scattering  sources  are  obviously  related  to  the  geometries  involved.  The 
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scattering  is  multiple  in  the  sense  that  the  various  wavelets  in  the  source  wavetrain 
generate  specific  Rg-wavelets.  Since  the  P-wave  energy  decrease  with  time  we  have  a 
corresponding  weakening  of  the  Rg-pulses. 

Further  into  the  coda,  say  30-70  s,  pulses  with  Sg  and  Sn  velocities  become  dominant. 
The  area  of  potential  scattering  sources  increases  greatly,  but  number  of  significant 
scatters  remain  modest.  Corresponding  bearings  are  mainly  from  west  to  northeast,  that  is 
the  most  mountainous  part  of  S.  Norway,  but  also  in  the  event  azimuth  direction.  The  S- 
wavelets  are  most  probably  generated  by  P-to-S  scattering  so  their  later  arrivals  simply 
reflect  longer  travelling  distances. 

In  our  analysis  we  have  found  very  little  evidence  of  P-to-P  scattering  (6-8  km/s).  We 
take  this  to  imply  that  structural  heterogeneities  in  the  crust/lithosphere  are  not  prominent 
enough  to  act  as  efficient  scattering  sources.  Even  the  nearby  Oslo  Rift  with  an  elevated 
Moho  of  3-5  km  (Kinck  et  al,  1990)  seemingly  do  not  generate  significant  P-to-P 
scattering.  Weak  lithospheric  P-scattering  at  depths  around  100  km  as  reported  by 
Troitskiy  et  al  (1981)  for  NORSAR  (100  km  aperture)  are  not  easily  separated  from 
source-end  contributions  for  NORESS  (only  3  km  aperture).  Also,  why  P-to-Rg,  and  P-to- 
S  conversions  are  far  more  efficient  than  P-to-P  conversions  remains  somewhat  puzzling 
but  these  problems  would  be  explored  using  3D  wavefield  synthesis. 

The  above  results  and  the  suggested  mechanism  for  coda  generation  strictly  apply  only  to 
NORESS  recordings  at  teleseismic  ranges.  It  was  rather  surprising  to  find  that  about  SO¬ 
SO  %  of  the  coda  (Table  2)  is  tied  to  wavelets  with  teleseismic  velocities  (the  coherent 
part)  while  10-30  %  appears  to  be  scattered  by  rough  topography  in  the  two  mentioned 
areas  near  the  array.  The  remaining  part  is  considered  to  be  diffuse  scattering  at  the 
receiver-end.  For  local  and  regional  distances  other  scattering  mechanisms  may  be  more 
relevant  since  we  would  be  dealing  with  horizontal  travelling  waves  of  shorter 
wavelengths.  An  other  interesting  observation  is  that  Rg-  propagation  seems  to  be  far 
more  efficient,  up  to  300  km,  in  the  vicinity  of  arrays  in  N.  Norway  (ARCESS)  and 
Finland  (FINES A)  according  to  Toks0z  et  al  (1990). 


5.2  Conclusions 
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In  this  study  we  have  examined  teleseismic  P-coda  waves  as  recorded  by  the  NORESS 
array  in  SE  Norway  in  the  time  interval  75  s  after  P-onset.  The  essential  elements  in  the 
data  analysis  are  as  follow: 

In  3C  analysis  the  phase  shift  attribute  (Z/R)  has  proved  efficient  in  indicating 
where  interference  takes  place. 

Semblance  analysis  and  beamforming  proved  very  valuable  in  estimating  the 
slowness  vector  from  non-P  phases  like  Rg,  Sg  and  Sn. 

Ideally,  3C  NORESS  recordings  would  be  more  useful  for  coda  studies  given 
additional  3C  stations. 

The  main  results  obtained  are  as  follow: 

Early  P-coda  at  NORESS  shows  prominent  P-to-Rg  scattering  from  local  hills  with 
pronounced  topographic  relief,  namely  Bronkeberget  (Dist~10  km,  Azi~80°)  and 
Skreikampen  (Dist~30  km,  Azi-ZZS^*).  Rg-scattering  is  both  forward  and  backward. 
P-to-S  scattering  is  also  quite  efficient,  continues  further  into  the  coda  than  Rg, 
but  scattering  locations  are  more  obscure. 

P-to-P  scattering  takes  place  in  the  immediate  vicinity  of  the  array  (within  4-5  s 
from  onset)  particularly  to  the  east,  but  is  generally  weakly  developed. 

Source-end  scattering  could  not  be  separated  from  extended  source  pulse  duration 
per  se. 

Scattering  is  multiple  in  the  sense  that  both  the  first  and  later  P-phases  from  the 
source  region  result  in  Rg-scattering  mainly  from  the  two  mentioned  locations. 

This  probably  applies  also  to  P-to-S  scattering. 
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Table  captions. 


Table  1:  Hypocenter  parameters  of  the  events  analyzed  (from  ISC  and  PDE  catalogues); 
distances  ard  azimuths  (in  degrees)  are  calculated  relative  to  the  NORESS  array.  Event 
regions  are:  KAM  -  Kamchatka,  JAP  -  Japan  (south  of  Honshu),  SEM  -  Semipalatinsk 
(East  Kazakh),  EUCl  and  HK2  -  Hindu  Kush,  MAD  -  NW  of  Madagascar,  LEE  - 
Leeward  Island,  ALA  -  Alaska. 


Table  2:  Shown  numbers  are  percentages  of  coherent  power  on  the  teleseismic  beam  for 
different  time  intervals.  Coherent  power  is  calculated  as  the  energy  on  the  beam  divided 
by  the  average  energy  of  the  single  channel  sensors  for  the  same  time  interval.  Data  used 
are  from  the  AO,  C-  and  D-ring  instruments  and  were  filtered  in  the  1-3  Hz  passband. 
The  beam  steering  directions  were  those  giving  the  highest  coherency  in  the  first  3  s  of 
the  signal. 
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Figure  captions. 


Fig.  1.  The  general  NORESS  array  area  with  scattering  locations  hatched.  The  arrows 
gives  the  direction  of  approach  of  the  events  analyzed  (Table  1).  The  array  configuration 
is  given  in  the  lower  right  insert  where  encircled  stations  indicate  3-component  sites. 

Fig.  2a.  Results  from  3-component  analysis  of  the  first  4  events  in  Table  1;  Each  event 
is  shown  in  a  compartment  box  of  4  sections,  namely; 

i)  The  3C  seismic  traces,  rotated  into  a  vertical-radial-tangential  (Z,R,T)  coordinate 
system. 

ii)  ML  probabilities  exceeding  0.50  of  a  P-wave  polarization.  For  all  events  the 
sliding  window  was  1.6s  long,  incremented  in  0.2s  in  time  and  4°  in  azimuth. 

iii)  The  phase  difference  between  the  radial  (R)  and  vertical  (Z)  components.  TTie 
particle  motion  is  pure  P,  elliptical  or  pure  SV  for  phase  differences  of  180°,  ± 

.  90°  and  0°,  respectively.  For  cases  where  the  arriving  wavelets  are  not 
approximately  in  the  radial  plane,  the  phase  difference  estimate  would  be 
unreliable. 

iv)  Apparent  3C-velocities  including  error  bars  as  calculated  from  angle  of  incidence 
estimates. 

Fig.  2b.  Results  from  3-component  (3C)  analysis  of  the  last  4  events  in  Table  1. 
Otherwise  caption  as  in  Fig.  2a. 

Fig.  3.  RMS  amplitudes  (counts)  for  average  single  channel  and  beam  traces  for  the 
events  analyzed  (Table  1).  The  traces  from  the  AO,  C-  and  D-ring  sensors  were  bandpass 
filtered  1-3  Hz  prior  to  RMS  calculations  in  a  4  s  window.  The  upper  traces  are  the 
average  single  channel  amplitudes  while  the  beam  amplitudes  are  shown  with  shading. 
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Fig.  4.  Outcome  of  residual  trace  semblance  analysis  aimed  at  identifying  Rg-wavelets  in 
the  coda.  Semblance  is  shown  as  a  function  time  and  azimuth  for  a  fixed  phase  velocity 
of  3.0  km/s.  The  calculation  of  residual  traces  are  explained  in  the  text.  All  traces  were 
filtered  in  the  1-3  Hz  passband  and  semblance  calculations  are  for  a  1.5  s  window.  The 
semblance  contouring  levels  start  at  0.25  and  increase  in  steps  of  0.05.  The  thin 
horizontal  lines  mark  azimuths  of  80°  and  225°,  while  the  thin  vertical  lines  indicate  P- 
onset  time. 

Fig.  5.  Examples  of  waveforms  estimated  through  beamforming.  For  each  event  the  ’best’ 
teleseismic  beam  is  shown  first  followed  by  Rg  waveforms  estimated  by  residual  trace 
beamforming  aimed  at  the  two  scattering  areas.  All  traces  were  filter  in  the  1-3  Hz 
passband. 

Fig.  6.  Coda  semblance  analysis  aimed  at  identifying  Rg-wavelets  in  a  40  sec  time 
segment  after  that  in  Fig.  4.  The  thin  horizontal  lines  mark  azimuths  of  80°  and  225°, 
while  the  thin  vertical  lines  indicate  40  s  after  P-onset  time.  Caption  otherwise  as  in  Fig. 
4. 

Fig.  7.  Coda  semblance  analysis  aimed  at  identifying  Sn-wavelets  for  the  same  time 
segment  as  in  Fig.  6.  The  fixed  phase  velocity  is  here  4.5  km/s  -  caption  otherwise  as  in 
Fig.  6.  In  cases  where  semblance  contours  coincide  with  those  in  Fig.  6,  the  true  phase 
velocity  is  probably  of  the  intermediate  Sg-type,  that  is  ca  3.6  km/s. 
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Event  1 

Date 

Origin-time 

Lat . 

Long. 

Dist . 

Azi . 

Dep. 

Mb/Ms 

KAM 

11 

Jan 

88 

21:07:29.7 

54.78N 

161. 66E 

62.1 

19.0 

43 

5. 8/4. 9 

JAP 

12 

Dec 

87 

04:51:50.5 

29.69N 

140. 02E 

80.3 

43.6 

164 

6.3/  - 

SEM 

14 

Sep 

88 

04:00:00.0 

49.87N 

78.82E 

37.9 

75.2 

0 

6 . 0/4 . 5 

HKl 

29 

Jul 

85 

07:54:44.5 

36.21N 

70.90E 

44.2 

95.7 

102 

6.5/  - 

HK2 

7 

May 

86 

23:25:25.9 

36.37N 

70.71E 

44.0 

95.8 

223 

5.6/  - 

MAD 

14 

May 

85 

13:25:01.2 

10.72S 

41.26E 

75.2 

149.8 

37 

5 . 8/5 . 6 

LEE 

16 

Mar 

85 

14:54:43.7 

16.98N 

62.46W 

67.4 

264.6 

20 

6 . 1/6 . 4 

ALA 

6 

Mar 

88 

22:35:38.1 

56.95N 

143. 03W 

60.6 

344.4 

10 

6 . 8/7 . 6 

Table:  1 


Time  after  P  onset  (s): 


Event 

0-5 

5-10 

10-35 

35-75 

KAM 

94 

63 

44 

44 

JAP 

94 

67 

34 

21 

SEM 

96 

63 

45 

42 

HKl 

92 

90 

63 

43 

HK2 

92 

52 

26 

28 

MAD 

97 

84 

75 

26 

LEE 

98 

95 

90 

71 

ALA 

93 

94 

61 

45 

Table:  2 
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RESEARCH  IMPLICATIONS  -  FUTURE  WORK 

We  are  a  bit  excited  about  the  new  scattering  results  (Sec.4)  and  have  already  taken  the 
first  steps  towards  similar  studies  for  the  ARCESS  (N  Norway)  and  FINESA  (Finland) 
arrays.  A  far  more  difficult  study  now  being  considered,  namely  to  compute  the  response 
of  the  NORESS  scattering  ’hills’  on  the  P-wave  recordings  using  finite  element  techniques. 

The  NORESS  and  AkCESS  arrays  aie  equipped  with  feur  3C  stations,  and  we  plan  to 
undenake  a  development  of  techniques  for  joint  processing  of  multichannel  3C  records.  The 
rationale  here  is  to  extract  information  about  wavetype  being  dominant  on  the  horizontal 
components.  We  would  continue  the  work  on  3C  signal  detection  where  preliminary  results 
so  far  are  promising.  For  example,  for  a  time  span  of  16  hours  a  single  3C  station  within 
NORESS  (C4)  detected  24  local/regional  events  of  which  21  were  ’shared’  by  the  array 
itself.  The  3C  detector  now  being  tested,  operates  on  low  thresholds  so  a  relative  large  set 
of  rules  are  being  established  for  eliminating  false  detections  and  besides  associate 
secondary  coda-detections  into  event  families. 

Finally,  some  of  the  above  research  work  is  undertaken  jointly  with  visiting  scientists  and 
colleagues  in  Norway,  Sweden  and  Finland. 
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